Skip to content

HTTPS clone URL

Subversion checkout URL

You can clone with
or
.
Download ZIP
Browse files

Bounding box functionality to e3d_bv and a e3d_qbvh for OpenCL usage.

  • Loading branch information...
commit f0f6a07287f488d08a301171d7cebecfb42e2e5a 1 parent 523822b
@dgud dgud authored
View
1  e3d/Makefile
@@ -34,6 +34,7 @@ MODULES= \
e3d_util \
e3d_vec \
e3d_q \
+ e3d_qbvh \
e3d_bv \
e3d__tri_quad \
e3d__meshclean \
View
9 e3d/e3d.hrl
@@ -13,6 +13,7 @@
%% 3D vector or location.
-type e3d_vector() :: {float(),float(),float()}.
+-type e3d_point() :: {float(),float(),float()}.
%% Compact 4x4 matrix representation.
-type e3d_compact_matrix() ::
@@ -27,6 +28,13 @@
float(),float(),float(),float(),
float(),float(),float(),float(),
float(),float(),float(),float()}.
+
+
+%% Types for e3d_bv
+-define(E3D_INFINITY, 1.79769313e308).
+-type e3d_bbox() :: {e3d_point(), e3d_point()}.
+-type e3d_bsphere() :: {e3d_point(), number()}.
+-type e3d_bv() :: e3d_bbox() | e3d_bsphere().
-record(e3d_face,
{vs=[], %List of vertex indices.
@@ -63,3 +71,4 @@
dir %Directory for file.
}).
+
View
277 e3d/e3d_bv.erl
@@ -2,9 +2,9 @@
%% e3d_bv.erl --
%%
%% Bounding volume operations.
-%% Currently only quickhull, and eigen-vecs calculation implemented
+%% Bounding boxes and quickhull, and eigen-vecs calculation implemented
%%
-%% Copyright (c) 2001-2008 Dan Gudmundsson
+%% Copyright (c) 2001-2010 Dan Gudmundsson
%%
%% See the file "license.terms" for information on usage and redistribution
%% of this file, and for a DISCLAIMER OF ALL WARRANTIES.
@@ -13,13 +13,282 @@
%%
-module(e3d_bv).
--export([eigen_vecs/1,quickhull/1,covariance_matrix/1]).
--import(e3d_vec, [dot/2,add/2,average/1,normal/1]).
+%% Bounding Box/Sphere
+-export([box/0,box/1,box/2,box/3,
+ union/2,
+ center/1,max_extent/1,
+ surface_area/1,volume/1,
+ sphere/1,
+ inside/2]).
+
+%% Other Stuff
+-export([eigen_vecs/1,quickhull/1,covariance_matrix/1]).
+-import(e3d_vec, [dot/2,add/2,average/1,average/2,dist_sqr/2,normal/1]).
-import(lists, [foldl/3]).
+-include("e3d.hrl").
+
-compile(inline).
+-define(BB_MIN, {{?E3D_INFINITY,?E3D_INFINITY,?E3D_INFINITY},
+ {-?E3D_INFINITY,-?E3D_INFINITY,-?E3D_INFINITY}}).
+
+
+%%--------------------------------------------------------------------
+%% @doc Creates a bounding box
+%% Infinite if no arguments is given
+%% Enclosing the two points or list if given.
+%% The box is expanded with Epsilon in all directions if given.
+%% box() -> infinite_box;
+%% box(point, point |[Epsilon])
+%% box([points]|[Epsilon])
+%% @end
+%%--------------------------------------------------------------------
+
+-spec box() -> e3d_bbox().
+box() ->
+ ?BB_MIN.
+
+-spec box([e3d_point()]) -> e3d_bbox().
+box([{X,Y,Z}|Vs]) ->
+ bounding_box_1(Vs, X, X, Y, Y, Z, Z).
+
+-spec box(e3d_point()|[e3d_point()], e3d_point()|float()) -> e3d_bbox().
+box([{X,Y,Z}|Vs], Expand) ->
+ {Min, Max} = bounding_box_1(Vs, X, X, Y, Y, Z, Z),
+ {add(Min, {-Expand,-Expand,-Expand}), add(Max,{Expand,Expand,Expand})};
+
+box({V10,V11,V12}, {V20,V21,V22}) ->
+ {MinX, MaxX} = if V10 < V20 -> {V10,V20};
+ true -> {V20,V10}
+ end,
+ {MinY, MaxY} = if V11 < V21 -> {V11, V21};
+ true -> {V21, V11}
+ end,
+ {MinZ, MaxZ} =if V12 < V22 -> {V12, V22};
+ true -> {V22, V12}
+ end,
+ {{MinX,MinY,MinZ},{MaxX,MaxY,MaxZ}}.
+
+-spec box(e3d_vector(), e3d_vector(), float()) -> e3d_bbox().
+box({V10,V11,V12}, {V20,V21,V22}, Expand) ->
+
+ {MinX, MaxX} = if V10 < V20 -> {V10,V20};
+ true -> {V20,V10}
+ end,
+ {MinY, MaxY} = if V11 < V21 -> {V11, V21};
+ true -> {V21, V11}
+ end,
+ {MinZ, MaxZ} =if V12 < V22 -> {V12, V22};
+ true -> {V22, V12}
+ end,
+ {add({MinX,MinY,MinZ}, {-Expand,-Expand,-Expand}),
+ add({MaxX,MaxY,MaxZ}, {Expand,Expand,Expand})}.
+
+%%--------------------------------------------------------------------
+%% @doc Creates the union of a bounding box and point| bounding box
+%% @end
+%%--------------------------------------------------------------------
+
+-spec union(e3d_bbox(), e3d_vector() | e3d_bbox()) -> e3d_bbox().
+union(BBox1 = {Min1={V10,V11,V12}, Max1={V20,V21,V22}},
+ BBox2 = {Min2={V30,V31,V32}, Max2={V40,V41,V42}}) ->
+ %% Avoid tuple construction if unnecessary
+ %% {{erlang:min(V10,V30), erlang:min(V11,V31), erlang:min(V12,V32)},
+ %% {erlang:max(V20,V40), erlang:max(V21,V41), erlang:max(V22,V42)}};
+ %% Bjorn fix the compiler :-)
+ %% The compiler can not optimize away the tuple construction
+ %% that's why the code looks like this.
+ if V10 =< V30 ->
+ if V11 =< V31 ->
+ if V12 =< V32 ->
+ if V20 >= V40 ->
+ if V21 >= V41 ->
+ if V22 >= V42 -> BBox1;
+ true -> {Min1, {V20,V21,V42}}
+ end;
+ true -> {Min1, {V20, V41, erlang:max(V22,V42)}}
+ end;
+ true ->
+ {Min1, {V40, erlang:max(V21,V41), erlang:max(V22,V42)}}
+ end;
+ true ->
+ {{V10,V11,V32}, max_point(Max1, Max2)}
+ end;
+ true ->
+ {{V10, V31, erlang:min(V12,V32)}, max_point(Max1, Max2)}
+ end;
+ true ->
+ if V31 =< V11 ->
+ if V32 =< V12 ->
+ if V40 >= V20 ->
+ if V41 >= V21 ->
+ if V42 >= V22 -> BBox2;
+ true -> {Min2, {V40,V41,V22}}
+ end;
+ true ->
+ {Min2, {V40, V21, erlang:max(V42,V22)}}
+ end;
+ true ->
+ {Min2, {V40, erlang:max(V41,V21), erlang:max(V42,V22)}}
+ end;
+ true ->
+ {{V30, V31, V12}, max_point(Max1, Max2)}
+ end;
+ true ->
+ {{V30, V11, erlang:min(V12, V32)}, max_point(Max1, Max2)}
+ end
+ end;
+
+union(BBox = {Min0={V10,V11,V12}, Max0 = {V20,V21,V22}},
+ Point = {V30,V31,V32}) ->
+ if V10 =< V30 ->
+ if V11 =< V31 ->
+ if V12 =< V32 ->
+ if V20 >= V30 ->
+ if V21 >= V31 ->
+ if V22 >= V32 -> BBox;
+ true -> {Min0, {V20,V21,V32}}
+ end;
+ true -> {Min0, {V20, V31, erlang:max(V22,V32)}}
+ end;
+ true ->
+ {Min0, {V30, erlang:max(V21,V31), erlang:max(V22,V32)}}
+ end;
+ true ->
+ {{V10,V11,V32}, max_point(Max0, Point)}
+ end;
+ true ->
+ {{V10, V31, erlang:min(V12,V32)}, max_point(Max0, Point)}
+ end;
+ true ->
+ if V31 =< V11 ->
+ if V32 =< V12 -> {Point, max_point(Max0, Point)};
+ true -> {{V30,V31,V12}, max_point(Max0, Point)}
+ end;
+ true -> {{V30,V11,erlang:min(V12,V32)}, max_point(Max0, Point)}
+ end
+ end.
+
+%%--------------------------------------------------------------------
+%% @doc Creates a bounding sphere from a bounding box
+%% @end
+%%--------------------------------------------------------------------
+
+-spec sphere(e3d_bbox()) -> e3d_bsphere().
+sphere(BB = {{_,_,_}, Max = {_,_,_}}) ->
+ Center = center(BB),
+ {Center,
+ case inside(BB, Center) of
+ true -> dist_sqr(Center, Max);
+ false -> 0.0
+ end}.
+
+%%--------------------------------------------------------------------
+%% @doc Returns the center of the bounding volume
+%% @end
+%%--------------------------------------------------------------------
+-spec center(e3d_bv()) -> e3d_point().
+center({Min = {_,_,_}, Max = {_,_,_}}) ->
+ average(Min,Max);
+center({Center, DistSqr}) when is_list(DistSqr) ->
+ Center.
+
+%%--------------------------------------------------------------------
+%% @doc Returns the surface area of the bounding volume
+%% @end
+%%--------------------------------------------------------------------
+-spec surface_area(e3d_bv()) -> float().
+surface_area({Min = {Minx,_,_}, Max = {MaxX,_,_}}) ->
+ if Minx > MaxX -> 0;
+ true ->
+ {X,Y,Z} = e3d_vec:sub(Max, Min),
+ X*Y+Y*Z+Z*X*2
+ end.
+
+
+%%--------------------------------------------------------------------
+%% @doc Returns the volume of the bounding volume
+%% @end
+%%--------------------------------------------------------------------
+-spec volume(e3d_bv()) -> float().
+volume({Min = {Minx,_,_}, Max = {MaxX,_,_}}) ->
+ if Minx > MaxX -> 0;
+ true ->
+ {X,Y,Z} = e3d_vec:sub(Max, Min),
+ X*Y*Z
+ end.
+
+
+%%--------------------------------------------------------------------
+%% @doc Returns true if point is inside baounding volume
+%% @end
+%%--------------------------------------------------------------------
+-spec inside(e3d_bv(), e3d_vector()) -> boolean().
+inside({{V10,V11,V12}, {V20,V21,V22}}, {V30,V31,V32}) ->
+ V10 >= V30 andalso V30 >= V20 andalso
+ V11 >= V31 andalso V31 >= V21 andalso
+ V12 >= V32 andalso V32 >= V22;
+inside({Center, DistSqr}, Point) when is_number(DistSqr) ->
+ dist_sqr(Center, Point) =< DistSqr.
+
+%%--------------------------------------------------------------------
+%% @doc Returns the largest dimension of the bounding box,
+%% 1 = X, 2 = Y, Z = 3
+%% @end
+%%--------------------------------------------------------------------
+
+-spec max_extent(e3d_bbox()) -> 1 | 2 | 3.
+
+max_extent({Min, Max}) ->
+ {X,Y,Z} = e3d_vec:sub(Max, Min),
+ if X > Y, X > Z -> 1;
+ Y > Z -> 2;
+ true -> 3
+ end.
+
+
+
+
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%% Internals
+max_point(Max0 = {V20,V21,V22}, Max1 = {V30,V31,V32}) ->
+ if V20 >= V30 ->
+ if V21 >= V31 ->
+ if V22 >= V32 -> Max0;
+ true -> {V20,V21,V32}
+ end;
+ true -> {V20, V31, erlang:max(V22,V32)}
+ end;
+ true ->
+ if V31 >= V21 ->
+ if V32 >= V22 -> Max1;
+ true -> {V30,V31,V22}
+ end;
+ true -> {V30, V21, erlang:max(V22,V32)}
+ end
+ end.
+
+bounding_box_1([{X,_,_}|_]=Vs, X0, X1, Y0, Y1, Z0, Z1) when X < X0 ->
+ bounding_box_1(Vs, X, X1, Y0, Y1, Z0, Z1);
+bounding_box_1([{X,_,_}|_]=Vs, X0, X1, Y0, Y1, Z0, Z1) when X > X1 ->
+ bounding_box_1(Vs, X0, X, Y0, Y1, Z0, Z1);
+bounding_box_1([{_,Y,_}|_]=Vs, X0, X1, Y0, Y1, Z0, Z1) when Y < Y0 ->
+ bounding_box_1(Vs, X0, X1, Y, Y1, Z0, Z1);
+bounding_box_1([{_,Y,_}|_]=Vs, X0, X1, Y0, Y1, Z0, Z1) when Y > Y1 ->
+ bounding_box_1(Vs, X0, X1, Y0, Y, Z0, Z1);
+bounding_box_1([{_,_,Z}|_]=Vs, X0, X1, Y0, Y1, Z0, Z1) when Z < Z0 ->
+ bounding_box_1(Vs, X0, X1, Y0, Y1, Z, Z1);
+bounding_box_1([{_,_,Z}|_]=Vs, X0, X1, Y0, Y1, Z0, Z1) when Z > Z1 ->
+ bounding_box_1(Vs, X0, X1, Y0, Y1, Z0, Z);
+bounding_box_1([_|Vs], X0, X1, Y0, Y1, Z0, Z1) ->
+ bounding_box_1(Vs, X0, X1, Y0, Y1, Z0, Z1);
+bounding_box_1([], X0, X1, Y0, Y1, Z0, Z1) ->
+ {{X0,Y0,Z0},{X1,Y1,Z1}}.
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
eigen_vecs(Vs) ->
Fs = quickhull(Vs),
SymMat = covariance_matrix(Fs),
View
387 e3d/e3d_qbvh.erl
@@ -0,0 +1,387 @@
+%%
+%% e3d_qbvh.erl --
+%%
+%% Implements a quad-bvh. Quad Bounding volume hierarchies
+%% http://www.uni-ulm.de/fileadmin/website_uni_ulm/iui.inst.100/institut/Papers/QBVH.pdf
+%%
+%% This code is not intended to be fast in erlang (if intersection
+%% is ever implemented), it is intended to be fast when sent down
+%% OpenCL.
+%%
+%% The code is inspired from LuxRays variant.
+%%
+%% Copyright (c) 2010, Dan Gudmundsson
+%%
+%% See the file "license.terms" for information on usage and redistribution
+%% of this file, and for a DISCLAIMER OF ALL WARRANTIES.
+%%
+
+-module(e3d_qbvh).
+
+-export([init/1, init/2]).
+
+-export([ray/4, print_tree/4]). %% For testing purposes
+
+-include("e3d.hrl").
+
+-import(lists, [reverse/1]).
+
+-define(EPSILON, 0.00001).
+-define(MAX_PRIMS_PER_LEAF, 4).
+-define(FULL_SWEEP, 4*?MAX_PRIMS_PER_LEAF).
+-define(SKIP_FACTOR, 3).
+-define(NB_BINS, 8).
+
+-define(F32, 32/float-native).
+-define(I32, 32/native).
+
+-define(DEBUG, 1).
+
+-record(qnode,
+ {children, % 4 Children
+ % Child = Leaf -> <<IsLeaf:1, NoOfLeafs:4, StartLeaf:27>>
+ % Node -> child index
+ bboxes % 4 Bounding boxes
+ }).
+
+-record(qray,
+ {o, % Origin
+ d, % Direction vector
+ min, % MinT
+ max}). % MaxT
+-type qray() :: #qray{}.
+
+%%--------------------------------------------------------------------
+%% @doc Builds the qbvh tree
+%% Indata is a list of {NumberOfFaces, GetVs(FaceIndex)} tuples.
+%% GetVs should return {VertexPos1,VertexPos2,VertexPos2}
+%% Face indecies are 0 numbered.
+%% @end
+%%--------------------------------------------------------------------
+
+-spec init([{NoFs :: integer(), GetVs :: function()}]) ->
+ {BB::e3d_bbox(), QTree :: binary(), QTris :: binary(), Map :: term()}.
+init(Data) ->
+ init(Data, []).
+
+init(FaceData, Opts) ->
+ Eps = proplists:get_value(epsilon, Opts, ?EPSILON),
+ {PrimBBs, PrimCs, WorldBB, CentriodBB} = init_bb(FaceData, Eps),
+ NPrims = array:size(PrimBBs),
+
+ %% Init index map and add 3 primitives at the end for the last quad
+ PrimIndx0 = array:from_list(lists:seq(0, NPrims-1) ++ [0,0,0]),
+ Nodes0 = {array:new([{default, create_empty_node()}]), 0},
+ {PrimIndx, {Nodes,_NodeLen}} =
+ build_tree(0, NPrims, Nodes0, PrimIndx0,
+ PrimBBs, PrimCs, WorldBB, CentriodBB, -1, 0, 0),
+ %% io:format("~nNodes ~p: ~p~n",[NPrims, array:to_list(PrimIndx)]),
+ %% io:format("Tree ~p: ~p~n",[_N, array:to_list(Nodes)]),
+ %% print_tree(0, Nodes, PrimIndx, 1),
+
+ {Tree, {_N, Tris}} = swap_prims(Nodes, PrimIndx, FaceData),
+ {WorldBB, array:foldl(fun binary_node/3, <<>>, Tree), Tris, []}.
+
+%%--------------------------------------------------------------------
+%% @doc Creates a ray
+%% @end
+%%--------------------------------------------------------------------
+-spec ray(e3d_point(), e3d_vector(), float(), float()) -> qray().
+ray(Orig, Vector, MinT, MaxT) ->
+ #qray{o=Orig, d=Vector, min=MinT, max=MaxT}.
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+print_tree(Index, Nodes, PrimIndx, Level) ->
+ #qnode{children=Children, bboxes=BBs} = array:get(Index, Nodes),
+ [print_child(C, Nodes, PrimIndx, Level)
+ || C <- lists:zip(tuple_to_list(Children), tuple_to_list(BBs))].
+
+print_child({-1,_BB}, _, _, Level) ->
+ io:format("~*.1c empty~n", [Level,$:]);
+print_child({{true, _Quads, Start},BB}, _, PrimIndx, Level) ->
+ io:format("~*.1c ~p ~p ~p ~p ~s~n",
+ [Level, $:|[array:get(Start+I,PrimIndx) || I <- lists:seq(0,3)]] ++ [f(BB)]);
+print_child({Node,BB}, Nodes, PrimIndx, Level) ->
+ io:format("~*.1c ~p ~s~n", [Level,$:, Node, f(BB)]),
+ print_tree(Node, Nodes, PrimIndx, Level+2).
+
+f({B1,B2}) ->
+ "<" ++ f(B1) ++ "," ++ f(B2) ++ ">";
+f({X,Y,Z}) ->
+ io_lib:format("{~.2g, ~.2g, ~.2g}", [X,Y,Z]);
+f(X) when is_float(X) ->
+ io_lib:format("~.2g", [X]);
+f(X) when is_integer(X) ->
+ io_lib:format("~p", [X]).
+
+init_bb(Data, Eps) when is_list(Data) ->
+ init_bb(Data, Eps, [], [], e3d_bv:box(), e3d_bv:box()).
+
+init_bb([{NrFs, Fun}|Rest], Eps, BBs, Cs, WBB0, CBB0) ->
+ {PrimBBs, PrimCs, WorldBB, CentriodBB} =
+ init_bb_1(NrFs-1, Fun, Eps, BBs, Cs, WBB0, CBB0),
+ init_bb(Rest, Eps, PrimBBs, PrimCs, WorldBB, CentriodBB);
+init_bb([], _Eps, PrimBBs, PrimCs, WorldBB, CentriodBB) ->
+ {array:from_list(reverse(PrimBBs)), array:from_list(reverse(PrimCs)),
+ WorldBB, CentriodBB}.
+
+init_bb_1(FaceIndex, GetVs, Eps, BBs, Cs, WBB0, CBB0)
+ when FaceIndex >= 0 ->
+ {V1,V2,V3} = GetVs(FaceIndex),
+ BB = e3d_bv:box([V1,V2,V3], Eps),
+ Center = e3d_bv:center(BB),
+ WBB = e3d_bv:union(WBB0, BB),
+ CBB = e3d_bv:union(CBB0, Center),
+ init_bb_1(FaceIndex-1, GetVs, Eps, [BB|BBs], [Center|Cs], WBB, CBB);
+init_bb_1(_, _, _Eps, BBs, Cs, WBB0, CBB0) ->
+ {BBs, Cs, WBB0, CBB0}.
+
+build_tree(Start, End, Nodes, PrimIndx, _PrimBBs, _PrimCs, NodeBox,
+ _CentriodBB, Parent, Child, _Depth)
+ when (End - Start) =< ?MAX_PRIMS_PER_LEAF ->
+ {PrimIndx, create_tmp_leaf(Parent, Child, Start, End, NodeBox, Nodes)};
+build_tree(Start, End, Nodes0, PrimIndx0, PrimBBs, PrimCs,
+ NodeBox, CentriodBB = {Cmin,Cmax}, Parent, Child, Depth) ->
+ Step = if (End - Start) < ?FULL_SWEEP -> 1; true -> ?SKIP_FACTOR end,
+ %% Find axis to do the split
+ Axis = e3d_bv:max_extent(CentriodBB),
+ K0 = CminAxis = element(Axis, Cmin),
+ CdistAxis = (element(Axis, Cmax) - K0),
+ try ?NB_BINS / CdistAxis of
+ K1 ->
+ case Depth rem 2 of
+ 0 ->
+ {Curr, Nodes1} =
+ create_node(Parent, Child, NodeBox, Nodes0),
+ LeftChild = 0,
+ RightChild = 2;
+ 1 ->
+ Nodes1 = Nodes0,
+ Curr = Parent,
+ LeftChild = Child,
+ RightChild = Child+1
+ end,
+
+ MinBin = binning(Start,End,Step,Axis,K0,K1,
+ PrimIndx0,PrimCs,PrimBBs,
+ erlang:make_tuple(?NB_BINS, 0),
+ erlang:make_tuple(?NB_BINS, e3d_bv:box())),
+ %% The split plane coordinate is the coordinate of the end of
+ %% the chosen bin along the split axis
+ SplitPos = CminAxis + (MinBin + 1) * CdistAxis/?NB_BINS,
+ Empty = e3d_bv:box(),
+ {Store, PrimIndx1, LBB, LCBB, RBB, RCBB} =
+ split(Start, End, Start, SplitPos, Axis,
+ PrimIndx0, PrimBBs, PrimCs,
+ Empty, Empty, Empty, Empty),
+ {PrimIndx, Nodes} =
+ build_tree(Start, Store, Nodes1, PrimIndx1, PrimBBs, PrimCs,
+ LBB, LCBB, Curr, LeftChild, Depth+1),
+ build_tree(Store, End, Nodes, PrimIndx, PrimBBs, PrimCs,
+ RBB, RCBB, Curr, RightChild, Depth+1)
+ catch
+ _:_ ->
+ %% If the bbox is a point, create a leaf, hoping there are
+ %% not more than 64 primitives that share the same center.
+ case (End - Start) > 64 of
+ true -> exit({geometry_error,
+ too_many_faces_with_same_center});
+ false ->
+ {PrimIndx0, create_tmp_leaf(Parent, Child, Start, End,
+ NodeBox, Nodes0)}
+ end
+ end.
+
+create_node(Parent, _Child, _NodeBox, {Nodes, NNodes}) when Parent < 0 ->
+ {NNodes, {Nodes, NNodes+1}};
+create_node(Parent, Child, NodeBox, {Nodes, NNodes}) ->
+ #qnode{bboxes=BBs,children=Ch} = array:get(Parent, Nodes),
+ Updated = #qnode{bboxes = setelement(Child+1, BBs, NodeBox),
+ children = setelement(Child+1, Ch, NNodes)},
+ {NNodes, {array:set(Parent, Updated, Nodes), NNodes+1}}.
+
+-ifdef(DEBUG).
+create_tmp_leaf(Parent, Child, Start, End, NodeBox, {Nodes, NNodes})
+ when Parent < 0 ->
+ %% The leaf is directly encoded in the intermediate node.
+ create_tmp_leaf(0, Child, Start, End, NodeBox, {Nodes, NNodes});
+create_tmp_leaf(Parent, Child, Start, End, NodeBox, {Nodes, NNodes}) ->
+ #qnode{bboxes=BBs,children=Ch} = array:get(Parent, Nodes),
+ Quads = ((End - Start)+3) div 4,
+ Updated =
+ #qnode{bboxes = setelement(Child+1, BBs, NodeBox),
+ children = setelement(Child+1, Ch, {true, Quads, Start})},
+ {array:set(Parent, Updated, Nodes), NNodes}.
+-else.
+create_tmp_leaf(Parent, Child, Start, End, NodeBox, {Nodes, NNodes})
+ when Parent < 0 ->
+ %% The leaf is directly encoded in the intermediate node.
+ create_tmp_leaf(0, Child, Start, End, NodeBox, {Nodes, NNodes});
+create_tmp_leaf(Parent, Child, Start, End, NodeBox, {Nodes, NNodes}) ->
+ #qnode{bboxes=BBs,children=Ch} = array:get(Parent, Nodes),
+ Quads = ((End - Start)+3) div 4,
+ ChildData = (1 bsl 31) bor (((Quads-1) band 16#F) bsl 27)
+ bor (Start band 16#7ffffff),
+ Updated = #qnode{bboxes = setelement(Child+1, BBs, NodeBox),
+ children = setelement(Child+1, Ch, ChildData)},
+ {array:set(Parent, Updated, Nodes), NNodes}.
+-endif().
+
+create_empty_node() ->
+ Empty = e3d_bv:box(),
+ EmptyLeaf = -1,
+ #qnode{bboxes = {Empty,Empty,Empty,Empty},
+ children = {EmptyLeaf, EmptyLeaf, EmptyLeaf, EmptyLeaf}}.
+
+swap_prims(Nodes, PrimIndx, Data) ->
+ swap_prims(0, array:get(0, Nodes), 1, Nodes, {0, <<>>}, PrimIndx, Data).
+
+swap_prims(Index, N0=#qnode{children=Children}, I, Nodes0,
+ Prims0, PrimIndx, Data)
+ when I < 5 ->
+ case element(I, Children) of
+ -1 ->
+ swap_prims(Index, N0, I+1, Nodes0, Prims0, PrimIndx, Data);
+ Leaf = {true, _Quads, _Start} ->
+ {NewLeaf, Prims} = swap_leaf(Leaf, Prims0, PrimIndx, Data),
+ Node = setelement(I, Children, NewLeaf),
+ swap_prims(Index, N0#qnode{children=Node}, I+1, Nodes0,
+ Prims, PrimIndx, Data);
+ NodeIndx ->
+ {Nodes,Prims} =
+ swap_prims(NodeIndx, array:get(NodeIndx, Nodes0), 1,
+ Nodes0, Prims0, PrimIndx, Data),
+ swap_prims(Index, N0, I+1, Nodes, Prims, PrimIndx, Data)
+ end;
+swap_prims(Index, Node, 5, Nodes, Prims, _PrimIndx, _Data) ->
+ {array:set(Index, Node, Nodes), Prims}.
+
+swap_leaf({true, Quads, Start}, {NQuads, Prims0}, PrimIndx, Data) ->
+ Leaf = 16#80000000 bor (((Quads-1) band 16#F) bsl 27)
+ bor (NQuads band 16#7ffffff),
+ Prims = swap_leafs(Quads, NQuads, Start, Prims0, PrimIndx, Data),
+ {Leaf, {NQuads+Quads, Prims}}.
+
+swap_leafs(NbQuads, Num, Offset, Prims, PrimIndx, Data)
+ when NbQuads > 0 ->
+ QT = quad_tri(Offset, PrimIndx, Data),
+ swap_leafs(NbQuads-1, Num+1, Offset+4,
+ <<Prims/binary, QT/binary>>, PrimIndx, Data);
+swap_leafs(_NbQuads, _Num, _Offset, Prims, _PrimIndx, _Data) ->
+ Prims.
+
+quad_tri(Offset, PrimIndx, Data) ->
+ Prim1 = array:get(Offset+0, PrimIndx),
+ {{V00x,V00y,V00z},{V01x,V01y,V01z},{V02x,V02y,V02z}} =
+ get_vs(Prim1, Data),
+ Prim2 = array:get(Offset+1, PrimIndx),
+ {{V10x,V10y,V10z},{V11x,V11y,V11z},{V12x,V12y,V12z}} =
+ get_vs(Prim2, Data),
+ Prim3 = array:get(Offset+2, PrimIndx),
+ {{V20x,V20y,V20z},{V21x,V21y,V21z},{V22x,V22y,V22z}} =
+ get_vs(Prim3, Data),
+ Prim4 = array:get(Offset+3, PrimIndx),
+ {{V30x,V30y,V30z},{V31x,V31y,V31z},{V32x,V32y,V32z}} =
+ get_vs(Prim4, Data),
+
+ <<V00x:?F32,V10x:?F32,V20x:?F32,V30x:?F32,
+ V00y:?F32,V10y:?F32,V20y:?F32,V30y:?F32,
+ V00z:?F32,V10z:?F32,V20z:?F32,V30z:?F32,
+ (V01x-V00x):?F32, (V11x-V10x):?F32, (V21x-V20x):?F32, (V31x-V30x):?F32,
+ (V01y-V00y):?F32, (V11y-V10y):?F32, (V21y-V20y):?F32, (V31y-V30y):?F32,
+ (V01z-V00z):?F32, (V11z-V10z):?F32, (V21z-V20z):?F32, (V31z-V30z):?F32,
+ (V02x-V00x):?F32, (V12x-V10x):?F32, (V22x-V20x):?F32, (V32x-V30x):?F32,
+ (V02y-V00y):?F32, (V12y-V10y):?F32, (V22y-V20y):?F32, (V32y-V30y):?F32,
+ (V02z-V00z):?F32, (V12z-V10z):?F32, (V22z-V20z):?F32, (V32z-V30z):?F32,
+ (Prim1):?I32, (Prim2):?I32, (Prim3):?I32, (Prim4):?I32
+ >>.
+
+get_vs(Prim0, Data) ->
+ {Prim, GetVs} = locate_data(Prim0, Data),
+ GetVs(Prim).
+
+locate_data(Prim, [{Id, Data}|_])
+ when Prim < Id ->
+ {Prim, Data};
+locate_data(Prim, [{Id,_}|Rest]) ->
+ locate_data(Prim-Id, Rest).
+
+binary_node(_, #qnode{children={C0,C1,C2,C3},
+ bboxes={{B0min,B0max},
+ {B1min,B1max},
+ {B2min,B2max},
+ {B3min,B3max}
+ }}, Bin) ->
+ {V00x,V00y,V00z} = B0min, {V01x,V01y,V01z} = B0max,
+ {V10x,V10y,V10z} = B1min, {V11x,V11y,V11z} = B1max,
+ {V20x,V20y,V20z} = B2min, {V21x,V21y,V21z} = B2max,
+ {V30x,V30y,V30z} = B3min, {V31x,V31y,V31z} = B3max,
+ <<Bin/binary,
+ V00x:?F32,V10x:?F32,V20x:?F32,V30x:?F32,
+ V00y:?F32,V10y:?F32,V20y:?F32,V30y:?F32,
+ V00z:?F32,V10z:?F32,V20z:?F32,V30z:?F32,
+ V01x:?F32,V11x:?F32,V21x:?F32,V31x:?F32,
+ V01y:?F32,V11y:?F32,V21y:?F32,V31y:?F32,
+ V01z:?F32,V11z:?F32,V21z:?F32,V31z:?F32,
+ C0:?I32,C1:?I32, C2:?I32, C3:?I32>>.
+
+%% Binning is relative to the centroids bbox and to the
+%% primitives' centroid.
+binning(Start, End, Step, Axis,K0,K1,PrimIndx,PrimCs, PrimBB, Bins, BinBBs)
+ when Start < End ->
+ Index = array:get(Start, PrimIndx),
+ Id = erlang:min(?NB_BINS, ceil(K1 * (element(Axis, array:get(Index, PrimCs)) - K0))),
+ binning(Start + Step, End, Step, Axis,K0,K1,PrimIndx,PrimCs,PrimBB,
+ setelement(Id, Bins, 1+element(Id, Bins)),
+ setelement(Id, BinBBs, e3d_bv:union(element(Id,BinBBs), array:get(Index,PrimBB))));
+binning(_Start, _End, _Step, _Axis,_K0,_K1,_PrimIndx,_PrimCs, _PrimBB, BinsT, BinBBsT) ->
+ %% Evaluate Split
+ Bins = tuple_to_list(BinsT),
+ BinBBs = tuple_to_list(BinBBsT),
+ LeftEstimates = count_binning(Bins, BinBBs, 0, e3d_bv:box(), []),
+ RightEstimates = count_binning(reverse(Bins), reverse(BinBBs), 0, e3d_bv:box(), []),
+ find_minbox(0, reverse(LeftEstimates), tl(RightEstimates), ?E3D_INFINITY, undefined).
+
+count_binning([Bin|Bins], [BinBB|BinBBs], CurrNb0, CurrBB0, Est) ->
+ CurrNb = Bin + CurrNb0,
+ CurrBB = e3d_bv:union(CurrBB0, BinBB),
+ Area = e3d_bv:surface_area(CurrBB),
+ count_binning(Bins, BinBBs, CurrNb, CurrBB, [Area*CurrNb|Est]);
+count_binning([], [], _CurrNb, _CurrBB, Volumes) ->
+ Volumes.
+
+find_minbox(I, [Left|Ls], [Right|Rs], MinCost, Bin) ->
+ case Left + Right of
+ Cost when Cost < MinCost ->
+ find_minbox(I+1, Ls, Rs, Cost, I);
+ _ ->
+ find_minbox(I+1, Ls, Rs, MinCost, Bin)
+ end;
+find_minbox(_, _, [], _MinCost, Bin) ->
+ Bin.
+
+ceil(X) when is_float(X) ->
+ round(0.5+X).
+
+%% Split
+split(I, End, Store, SplitPos, Axis, PrimIndx0, PrimBBs, PrimCs, LBB0, LCBB0, RBB0, RCBB0)
+ when I < End ->
+ Index = array:get(I, PrimIndx0),
+ case element(Axis, array:get(Index, PrimCs)) =< SplitPos of
+ true ->
+ PrimIndx1 = array:set(I, array:get(Store, PrimIndx0), PrimIndx0),
+ PrimIndx = array:set(Store, Index, PrimIndx1),
+
+ LBB = e3d_bv:union(LBB0, array:get(Index, PrimBBs)),
+ LCBB = e3d_bv:union(LCBB0, array:get(Index, PrimCs)),
+ split(I+1, End, Store+1, SplitPos, Axis, PrimIndx, PrimBBs, PrimCs, LBB, LCBB, RBB0, RCBB0);
+ false ->
+ RBB = e3d_bv:union(RBB0, array:get(Index, PrimBBs)),
+ RCBB = e3d_bv:union(RCBB0, array:get(Index, PrimCs)),
+ split(I+1, End, Store, SplitPos, Axis, PrimIndx0, PrimBBs, PrimCs, LBB0, LCBB0, RBB, RCBB)
+ end;
+split(_I, _End, Store, _SplitPos, _Axis, PrimIndx, _PrimBBs, _PrimCs, LBB, LCBB, RBB, RCBB) ->
+ {Store, PrimIndx, LBB, LCBB, RBB, RCBB}.
+%%
+
View
24 e3d/e3d_vec.erl
@@ -260,27 +260,11 @@ average({V10,V11,V12}, {V20,V21,V22}, {V30,V31,V32}, {V40,V41,V42})
L = 0.25,
{L*(V10+V20+V30+V40),L*(V11+V21+V31+V41),L*(V12+V22+V32+V42)}.
+
+%% Should be removed and calls should be changed to e3d_bv instead.
-spec bounding_box([e3d_vector()]) -> [e3d_vector()].
-
-bounding_box([{X,Y,Z}|Vs]) ->
- bounding_box_1(Vs, X, X, Y, Y, Z, Z).
-
-bounding_box_1([{X,_,_}|_]=Vs, X0, X1, Y0, Y1, Z0, Z1) when X < X0 ->
- bounding_box_1(Vs, X, X1, Y0, Y1, Z0, Z1);
-bounding_box_1([{X,_,_}|_]=Vs, X0, X1, Y0, Y1, Z0, Z1) when X > X1 ->
- bounding_box_1(Vs, X0, X, Y0, Y1, Z0, Z1);
-bounding_box_1([{_,Y,_}|_]=Vs, X0, X1, Y0, Y1, Z0, Z1) when Y < Y0 ->
- bounding_box_1(Vs, X0, X1, Y, Y1, Z0, Z1);
-bounding_box_1([{_,Y,_}|_]=Vs, X0, X1, Y0, Y1, Z0, Z1) when Y > Y1 ->
- bounding_box_1(Vs, X0, X1, Y0, Y, Z0, Z1);
-bounding_box_1([{_,_,Z}|_]=Vs, X0, X1, Y0, Y1, Z0, Z1) when Z < Z0 ->
- bounding_box_1(Vs, X0, X1, Y0, Y1, Z, Z1);
-bounding_box_1([{_,_,Z}|_]=Vs, X0, X1, Y0, Y1, Z0, Z1) when Z > Z1 ->
- bounding_box_1(Vs, X0, X1, Y0, Y1, Z0, Z);
-bounding_box_1([_|Vs], X0, X1, Y0, Y1, Z0, Z1) ->
- bounding_box_1(Vs, X0, X1, Y0, Y1, Z0, Z1);
-bounding_box_1([], X0, X1, Y0, Y1, Z0, Z1) ->
- [{X0,Y0,Z0},{X1,Y1,Z1}].
+bounding_box(List) when is_list(List) ->
+ tuple_to_list(e3d_bv:box(List)).
-spec degrees(e3d_vector(), e3d_vector()) -> float().
View
139 e3d/test/qbvh.erl
@@ -0,0 +1,139 @@
+%%%-------------------------------------------------------------------
+%%% @author Dan Gudmundsson <dgud@erlang.org>
+%%% @copyright (C) 2010, Dan Gudmundsson
+%%% @doc Some test cases for qbvh
+%%%
+%%% @end
+%%% Created : 9 Jun 2010 by Dan Gudmundsson <dgud@erlang.org>
+%%%-------------------------------------------------------------------
+-module(qbvh).
+
+-compile(export_all).
+-include_lib("e3d.hrl").
+-include_lib("cl/include/cl.hrl").
+
+-define(F32, 32/float-native).
+-define(I32, 32/native).
+
+-define(VS, {{ 0.4, 1.0, -0.4}, %1
+ { 0.4, 0.0, -0.4}, %2
+ {-0.4, 0.0, -0.4},
+ {-0.4, 1.0, -0.4}, %4
+ {-0.4, 1.0, 0.4},
+ { 0.4, 1.0, 0.4}, %6
+ { 0.4, 0.0, 0.4},
+ {-0.4, 0.0, 0.4}}).%8
+
+-define(FS, [[0,1,2,3],
+ [2,7,4,3],
+ [0,5,6,1],
+ [5,4,7,6],
+ [5,0,3,4],
+ [6,7,2,1]]).
+
+-define(RAYHIT_SZ, 16).
+
+go() -> start().
+
+start() ->
+ Fs = lists:append([ [[V1,V2,V3],[V1,V3,V4]] || [V1,V2,V3,V4] <- ?FS]),
+ TFs = list_to_tuple(Fs),
+ format_faces(0, Fs, ?VS),
+ Qbvh = e3d_qbvh:init([{tuple_size(TFs),
+ fun(Face) ->
+ [V1,V2,V3] = element(Face+1, TFs),
+ {element(V1+1,?VS), element(V2+1,?VS), element(V3+1, ?VS)}
+ end}]),
+ {Rays, Res} = init_cl(Qbvh),
+ check_rays(Rays, Res, array:from_list(Fs), ?VS),
+ ok.
+
+check_rays([_Miss|Ds], <<_:12/binary, 16#FFFFFFFF:32, Rest/binary>>, Fs, Vs) ->
+ check_rays(Ds, Rest, Fs, Vs);
+check_rays([Dir|Ds], <<T:?F32,B1:?F32,B2:?F32,Face:?I32, Rest/binary>>, Fs, Vs) ->
+ io:format("~s => ~s ~s ~s ~s~n",[f(Dir), f(T),f(B1),f(B2),f(Face)]),
+ FVs = [V1,V2,V3] = array:get(Face, Fs),
+ io:format(" ~p => [~s,~s,~s]~n",
+ [FVs, f(element(V1+1,Vs)),f(element(V2+1,Vs)),f(element(V3+1,Vs))]),
+ check_rays(Ds, Rest, Fs, Vs);
+check_rays([],<<>>, _, _) ->
+ ok.
+
+format_faces(I, [FVs=[V1,V2,V3]|Fs], Vs) ->
+ io:format("~p ~p => [~s,~s,~s]~n",
+ [I, FVs, f(element(V1+1,Vs)),f(element(V2+1,Vs)),f(element(V3+1,Vs))]),
+ format_faces(I+1, Fs,Vs);
+format_faces(_,[],_) -> ok.
+
+
+f({X,Y,Z}) ->
+ io_lib:format("{~.2f, ~.2f, ~.2f}", [X,Y,Z]);
+f(X) when is_float(X) ->
+ io_lib:format("~.2f", [X]);
+f(X) when is_integer(X) ->
+ io_lib:format("~p", [X]).
+
+init_cl(Qbvh) ->
+ CL = clu:setup(all),
+ try
+ {ok,Name} = cl:get_platform_info(CL#cl.platform, name),
+ io:format("Platform: ~s~n~n", [Name]),
+ [Device] = CL#cl.devices,
+ {ok,DeviceInfo} = cl:get_device_info(Device),
+ ImageSupport = proplists:get_value(image_support, DeviceInfo, false),
+ MaxWorkGroupSize = proplists:get_value(max_work_group_size, DeviceInfo, 4),
+
+ %% io:format("DeviceInfo: ~p\n", [DeviceInfo])
+ %%Kernel = compile(CL, "test_kernel.cl"),
+ Kernel = compile(CL, "qbvh_kernel.cl"),
+ {ok,Queue} = cl:create_queue(CL#cl.context,Device,[]),
+ {ok,Local} = cl:get_kernel_workgroup_info(Kernel, Device, work_group_size),
+ io:format("work_group_size = ~p ~p\n", [Local, MaxWorkGroupSize]),
+
+ standard(CL#cl.context, Kernel, Queue, 64, Qbvh)
+ catch
+ _Error:Reason ->
+ io:format("Error ~p~n", [Reason])
+ after
+ io:format("CL ~p~n",[CL]),
+ clu:teardown(CL)
+ end.
+
+compile(CL, File) ->
+ {ok, Bin} = file:read_file(File),
+ {ok, Program} = clu:build_source(CL, Bin),
+ {ok, [Kernel]} = cl:create_kernels_in_program(Program),
+ Kernel.
+
+standard(Context, Kernel, Queue, WorkGroupSz, {_BB, QnodesBin, QtrisBin, _}) ->
+ {ok, QNodes} = cl:create_buffer(Context, [read_only, copy_host_ptr], byte_size(QnodesBin), QnodesBin),
+ {ok, QTris} = cl:create_buffer(Context, [read_only, copy_host_ptr], byte_size(QtrisBin), QtrisBin),
+
+ {RaysBin, Dirs, NoRays} = create_rays(),
+ {ok, Rays} = cl:create_buffer(Context, [read_only], byte_size(RaysBin), RaysBin),
+ {ok, Hits} = cl:create_buffer(Context, [write_only], NoRays * ?RAYHIT_SZ),
+
+ clu:apply_kernel_args(Kernel, [Rays,Hits,QNodes,QTris,NoRays, {local,24*WorkGroupSz*4}]),
+ %% {ok,Event} = cl:enqueue_task(Queue, Kernel, []),
+ {ok, Run} = cl:enqueue_nd_range_kernel(Queue, Kernel,
+ [NoRays], [WorkGroupSz], []),
+ {ok, Done} = cl:enqueue_read_buffer(Queue,Hits,0, NoRays*?RAYHIT_SZ, [Run]),
+ {ok, Result} = cl:wait(Done,1000),
+
+ {Dirs, Result}.
+
+create_rays() ->
+ Origo = {0.0, 0.5, 5.0},
+ MinT = 0.000005,
+ MaxT = ?E3D_INFINITY,
+ Dirs = [calc_dir(Dir/255) || Dir <- lists:seq(0, 255)],
+ {<< << (create_ray(Origo, Dir, MinT, MaxT))/binary >>
+ || Dir <- Dirs >>, Dirs, 256}.
+
+calc_dir(Dir) ->
+ e3d_vec:norm(e3d_q:vec_rotate({1.0,0.0,0.0}, e3d_q:from_angle_axis(Dir*360, {0.0,1.0,0.0}))).
+
+create_ray({OX,OY,OZ}, {DX,DY,DZ}, MinT, MaxT) ->
+ <<OX:?F32, OY:?F32, OZ:?F32,
+ DX:?F32, DY:?F32, DZ:?F32,
+ MinT:?F32, MaxT:?F32>>.
Please sign in to comment.
Something went wrong with that request. Please try again.