Skip to content

HTTPS clone URL

Subversion checkout URL

You can clone with HTTPS or Subversion.

Download ZIP
Browse files

Implemented intersection function in erlang. And fixed bugs.

I had to implement the intersection testing (called ray_trace
currently) in erlang to debug the strange renderings I got.

e3d_transform:lookat/3 fixed to behave with opengl standard matrices.

Changed ?E3D_INFINTY to be 32 bits float infinity, otherwise it can't
be moved to/from binaries. This also helps because you can safely
multiple values below that and don't crash because the float get's to
large.

Fixed a serious bug in bbox:union/2 code which was the main culprit
to my problems.
  • Loading branch information...
commit 8be7f1c7fc7356f6cd58f316598ffbf3b4097c99 1 parent 2665625
@dgud dgud authored
View
2  e3d/e3d.hrl
@@ -31,7 +31,7 @@
%% Types for e3d_bv
--define(E3D_INFINITY, 1.79769313e308).
+-define(E3D_INFINITY, 3.402823e+38). %% 32 bits float max
-type e3d_bbox() :: {e3d_point(), e3d_point()}.
-type e3d_bsphere() :: {e3d_point(), number()}.
-type e3d_bv() :: e3d_bbox() | e3d_bsphere().
View
2  e3d/e3d_bv.erl
@@ -131,7 +131,7 @@ union(BBox1 = {Min1={V10,V11,V12}, Max1={V20,V21,V22}},
{Min2, {V40, V21, erlang:max(V42,V22)}}
end;
true ->
- {Min2, {V40, erlang:max(V41,V21), erlang:max(V42,V22)}}
+ {Min2, {V20, erlang:max(V41,V21), erlang:max(V42,V22)}}
end;
true ->
{{V30, V31, V12}, max_point(Max1, Max2)}
View
433 e3d/e3d_qbvh.erl
@@ -4,9 +4,8 @@
%% 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.
+%% This code is not intended to be fast in erlang,
+%% it is intended to be fast when sent down OpenCL.
%%
%% The code is inspired from LuxRays variant.
%%
@@ -20,20 +19,22 @@
-export([init/1, init/2]).
--export([ray/4, print_tree/4]). %% For testing purposes
+-export([ray/2, ray/4, ray_trace/2, print_tree/3]). %% For testing purposes
-include("e3d.hrl").
-import(lists, [reverse/1]).
--define(EPSILON, 0.00001).
+-define(EPSILON, 0.0001).
-define(MAX_PRIMS_PER_LEAF, 4).
--define(FULL_SWEEP, 4*?MAX_PRIMS_PER_LEAF).
+-define(FULL_SWEEP, 8*?MAX_PRIMS_PER_LEAF).
-define(SKIP_FACTOR, 3).
-define(NB_BINS, 8).
-define(F32, 32/float-native).
--define(I32, 32/native).
+-define(I32, 32/signed-native).
+-define(QNODE_SZ, ((24+4)*4)).
+-define(QTRI_SZ, ((12*3+4)*4)).
-define(DEBUG, 1).
@@ -44,12 +45,12 @@
bboxes % 4 Bounding boxes
}).
--record(qray,
- {o, % Origin
- d, % Direction vector
- min, % MinT
- max}). % MaxT
--type qray() :: #qray{}.
+-record(ray,
+ {o,d, % Origin, Direction vector
+ n, f}). % Near far (or MinT MaxT)
+
+-record(hit, {t, b1, b2, f = 16#ffffffff}).
+-record(qtri, {v,e1,e2,f}).
%%--------------------------------------------------------------------
%% @doc Builds the qbvh tree
@@ -70,14 +71,18 @@ init(FaceData, Opts) ->
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]),
+ PIndx0 = array:from_list(lists:seq(0, NPrims-1)),
+ PIndx1 = array:set(NPrims, 0, PIndx0),
+ PIndx2 = array:set(NPrims+1, 0, PIndx1),
+ PIndx = array:set(NPrims+2, 0, PIndx2),
Nodes0 = {array:new([{default, create_empty_node()}]), 0},
{PrimIndx, {Nodes,_NodeLen}} =
- build_tree(0, NPrims, Nodes0, PrimIndx0,
+ build_tree(0, NPrims, Nodes0, PIndx,
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),
+
+%% io:format("~nNodes ~p: ~p~n",[NPrims, array:to_list(PrimIndx)]),
+%% io:format("Tree ~p: ~p~n",[_NodeLen, array:to_list(Nodes)]),
+%% print_tree(Nodes, PrimIndx, WorldBB),
{Tree, {_N, Tris}} = swap_prims(Nodes, PrimIndx, FaceData),
{WorldBB, array:foldl(fun binary_node/3, <<>>, Tree), Tris, []}.
@@ -86,44 +91,28 @@ init(FaceData, Opts) ->
%% @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}.
-%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+ray(Orig, Vector) ->
+ ray(Orig, Vector, ?EPSILON*10, ?E3D_INFINITY).
+ray(Orig, Vector, MinT, MaxT) ->
+ #ray{o=Orig, d=Vector, n=MinT, f=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).
+ray_trace(Ray = #ray{d=Dir}, {_BB, Tree, Tris, _Map}) ->
+%% io:format("Dir ~p ~n", [Ray]),
+ Swap = inv_sign(Dir),
+ ray_trace(Ray, [0], #hit{}, Swap, Tree, Tris).
-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(lists:reverse(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)),
+ {array:from_list(PrimBBs), array:from_list(PrimCs),
WorldBB, CentriodBB}.
init_bb_1(FaceIndex, GetVs, Eps, BBs, Cs, WBB0, CBB0)
@@ -162,7 +151,6 @@ build_tree(Start, End, Nodes0, PrimIndx0, PrimBBs, PrimCs,
LeftChild = Child,
RightChild = Child+1
end,
-
MinBin = binning(Start,End,Step,Axis,K0,K1,
PrimIndx0,PrimCs,PrimBBs,
erlang:make_tuple(?NB_BINS, 0),
@@ -175,6 +163,8 @@ build_tree(Start, End, Nodes0, PrimIndx0, PrimBBs, PrimCs,
split(Start, End, Start, SplitPos, Axis,
PrimIndx0, PrimBBs, PrimCs,
Empty, Empty, Empty, Empty),
+ %% io:format("~p ~f: ~n",[Axis, SplitPos]),
+ %% print_split(Start,Store,End,PrimIndx1),
{PrimIndx, Nodes} =
build_tree(Start, Store, Nodes1, PrimIndx1, PrimBBs, PrimCs,
LBB, LCBB, Curr, LeftChild, Depth+1),
@@ -201,7 +191,6 @@ create_node(Parent, Child, NodeBox, {Nodes, NNodes}) ->
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.
@@ -209,24 +198,9 @@ create_tmp_leaf(Parent, 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)},
+ children = setelement(Child+1, Ch, {true, Quads, Start})},
{array:set(Parent, Updated, Nodes), NNodes}.
--endif().
create_empty_node() ->
Empty = e3d_bv:box(),
@@ -284,17 +258,20 @@ quad_tri(Offset, PrimIndx, Data) ->
Prim4 = array:get(Offset+3, PrimIndx),
{{V30x,V30y,V30z},{V31x,V31y,V31z},{V32x,V32y,V32z}} =
get_vs(Prim4, Data),
-
+%% io:format("~w ",[[Prim1,Prim2,Prim3,Prim4]]),
<<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
+
+ (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) ->
@@ -331,17 +308,22 @@ binary_node(_, #qnode{children={C0,C1,C2,C3},
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))),
+ BinPos = ceil(K1 * (element(Axis, array:get(Index, PrimCs)) - K0)),
+ Id = erlang:min(?NB_BINS, BinPos),
+ BinBB = e3d_bv:union(element(Id,BinBBs), array:get(Index,PrimBB)),
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) ->
+ setelement(Id, BinBBs, BinBB));
+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).
+ 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,
@@ -365,7 +347,8 @@ ceil(X) when is_float(X) ->
round(0.5+X).
%% Split
-split(I, End, Store, SplitPos, Axis, PrimIndx0, PrimBBs, PrimCs, LBB0, LCBB0, RBB0, RCBB0)
+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
@@ -375,13 +358,301 @@ split(I, End, Store, SplitPos, Axis, PrimIndx0, PrimBBs, PrimCs, LBB0, LCBB0, RB
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);
+ 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)
+ 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) ->
+split(_I, _End, Store, _SplitPos, _Axis, PrimIndx, _PrimBBs, _PrimCs,
+ LBB, LCBB, RBB, RCBB) ->
{Store, PrimIndx, LBB, LCBB, RBB, RCBB}.
-%%
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%% DEBUGGING
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+is_leaf(Index) -> Index < 0.
+is_empty(Index) -> Index =:= 16#ffffffff.
+no_quads(Index) -> ((Index bsr 27) band 16#F) + 1.
+first_quad(Index) -> Index band 16#7ffffff.
+
+inv_sign({X,Y,Z}) ->
+ {inv(X),inv(Y),inv(Z),
+ X < 0.0, Y < 0.0, Z < 0.0}.
+inv(N) ->
+ try 1.0/N
+ catch _:_ ->
+ ?E3D_INFINITY
+ end.
+
+ray_trace(Ray, [NodeData|Stack], Hit, Swap, Tree, Tris) ->
+ case is_leaf(NodeData) of
+ false ->
+ Node = get_node(NodeData, Tree),
+ Visit = qbbox_intersect(Ray, Node, Swap),
+ S = add_stack(Visit, Node#qnode.children, Stack),
+ ray_trace(Ray, S, Hit, Swap, Tree, Tris);
+ true ->
+ QuadPrims = no_quads(NodeData),
+ Offset = first_quad(NodeData),
+ {R,H} = qtri_intersect(Offset, Offset+QuadPrims, Ray,Hit, Tris),
+ ray_trace(R,Stack,H,Swap,Tree,Tris)
+ end;
+ray_trace(_, [], Hit, _, _, _) -> Hit.
+
+qbbox_intersect(#ray{o={X,Y,Z}, n=N,f=F},
+ #qnode{bboxes={X0,Y0,Z0,X1,Y1,Z1}},
+ {InvX,InvY,InvZ,SwapX,SwapY,SwapZ}) ->
+ %% Swap signs
+ {Xmin,Xmax} = if SwapX -> {X1,X0}; true -> {X0,X1} end,
+ {Ymin,Ymax} = if SwapY -> {Y1,Y0}; true -> {Y0,Y1} end,
+ {Zmin,Zmax} = if SwapZ -> {Z1,Z0}; true -> {Z0,Z1} end,
+
+ Min0 = [N,N,N,N],
+ Max0 = [F,F,F,F],
+ Max = fun(In,C,MC,Inv) -> max(In, (MC-C)*Inv) end,
+ Min = fun(In,C,MC,Inv) -> min(In, (MC-C)*Inv) end,
+ %% X coord
+ MinX = apply2c2(Min0,X,Xmin,InvX,Max), % max(Min0, (Xmin-X)*InvX)
+ MaxX = apply2c2(Max0,X,Xmax,InvX,Min),
+ %%io:format("BBox X ~p~n", [apply2(MaxX, MinX, fun(MaxV,MinV) -> MaxV >= MinV end)]),
+ %% Y coord
+ MinY = apply2c2(MinX,Y,Ymin,InvY,Max),
+ MaxY = apply2c2(MaxX,Y,Ymax,InvY,Min),
+ %%io:format("BBox Y ~p~n", [apply2(MaxY, MinY, fun(MaxV,MinV) -> MaxV >= MinV end)]),
+ %% Z coord
+ MinZ = apply2c2(MinY,Z,Zmin,InvZ,Max),
+ MaxZ = apply2c2(MaxY,Z,Zmax,InvZ,Min),
+ %% io:format("BBox intersect ~p~n", [apply2(MaxZ, MinZ, fun(MaxV,MinV) -> MaxV >= MinV end)]),
+ apply2(MaxZ, MinZ, fun(MaxV,MinV) -> MaxV >= MinV end).
+
+add_stack([true|Visit], [Ch|Chs], Stack) ->
+ case is_empty(Ch) of
+ true -> add_stack(Visit,Chs,Stack);
+ false ->
+%% print_stack([Ch|Stack]),
+ add_stack(Visit,Chs,[Ch|Stack])
+ end;
+add_stack([_|Visit], [_|Chs], Stack) ->
+ add_stack(Visit,Chs,Stack);
+add_stack([],[],Stack) -> Stack.
+
+print_stack([Top|Rest]) ->
+ case is_leaf(Top) of
+ true ->
+ io:format(" leaf ~p(~p), ", [first_quad(Top), no_quads(Top)]);
+ false ->
+ io:format(" node ~p, ",[Top])
+ end,
+ print_stack(Rest);
+print_stack([]) ->
+ io:format("~n",[]).
+
+qtri_intersect(Offset, Size, Ray0, Hit0, Tris)
+ when Offset < Size ->
+ Tri = get_qtriangles(Offset, Tris),
+ {Ray,Hit} = qtri_intersect(Ray0, Tri, Hit0),
+ qtri_intersect(Offset+1, Size, Ray, Hit, Tris);
+qtri_intersect(_, _, Ray, Hit, _) -> {Ray, Hit}.
+
+qtri_intersect(#ray{d={DirX,DirY,DirZ},
+ o={Ox,Oy,Oz},
+ n=Near,f=Far} = Ray0,
+ #qtri{v = {Tx,Ty,Tz},
+ e1 = {E1X,E1Y,E1Z},
+ e2 = {E2X,E2Y,E2Z},
+ f = Fs},
+ Hit0) ->
+ %% Calc B1
+ Bary = fun(E1,D1,E2,D2) -> (D1*E1) - (D2*E2) end,
+ MulAdd = fun(SX,EX,SY,EY,SZ,EZ) ->
+ SX*EX+SY*EY+SZ*EZ
+ end,
+ MulAddDiv = fun(SX,EX,SY,EY,SZ,EZ,Div) ->
+ try (SX*EX+SY*EY+SZ*EZ)/Div
+ catch _:_ -> -1.0 end
+ end,
+
+ S1X = apply2c2(E2Z,DirY,E2Y,DirZ,Bary),
+ S1Y = apply2c2(E2X,DirZ,E2Z,DirX,Bary),
+ S1Z = apply2c2(E2Y,DirX,E2X,DirY,Bary),
+
+ Div = apply6(S1X,E1X,S1Y,E1Y,S1Z,E1Z,MulAdd),
+
+ Dx = [Ox-X || X <- Tx],
+ Dy = [Oy-Y || Y <- Ty],
+ Dz = [Oz-Z || Z <- Tz],
+
+ B1 = apply7(Dx,S1X,Dy,S1Y,Dz,S1Z,Div,MulAddDiv),
+ %% Calc B2
+
+ S2X = apply4(E1Z,Dy,E1Y,Dz,Bary),
+ S2Y = apply4(E1X,Dz,E1Z,Dx,Bary),
+ S2Z = apply4(E1Y,Dx,E1X,Dy,Bary),
+ B2 = apply4c3(S2X,DirX,S2Y,DirY,S2Z,DirZ,Div,MulAddDiv),
+
+ %% Calc. B0
+ B0 = apply2(B1, B2, fun(Bx,By) -> 1.0 - Bx - By end),
+
+ %% Calc T
+ T = apply7(E2X,S2X,E2Y,S2Y,E2Z,S2Z,Div, MulAddDiv),
+
+ case check_qtris(Div,B0,B1,B2,T,Fs,Near,Far,undefined) of
+ undefined ->
+ {Ray0,Hit0};
+ Hit ->
+% io:format(" Hit ~p~n", [Hit]),
+ {Ray0#ray{f=Hit#hit.t},Hit}
+ end.
+
+check_qtris([D|Ds],[B0|B0s],[B1|B1s],[B2|B2s],[T|Ts],[Id|Ids],N,F,Prev) ->
+ if D =/= 0.0, B0 >= 0.0, B1 >= 0.0, B2 >= 0.0,
+ T > N, T < F ->
+ Hit = #hit{t=T,b1=B1,b2=B2,f=Id},
+ check_qtris(Ds,B0s,B1s,B2s,Ts,Ids,N,T,Hit);
+ true ->
+ check_qtris(Ds,B0s,B1s,B2s,Ts,Ids,N,F,Prev)
+ end;
+check_qtris([],[],[],[],[],[],_,_,Hit) -> Hit.
+
+%% Helpers
+
+apply2([A,B,C,D],[X,Y,Z,W], Fun) ->
+ [Fun(A,X), Fun(B,Y), Fun(C,Z), Fun(D,W)].
+
+apply2c2([A,B,C,D],R, [X,Y,Z,W], I, Fun) ->
+ [Fun(A,R,X,I), Fun(B,R,Y,I), Fun(C,R,Z,I), Fun(D,R,W,I)].
+
+apply4([A1,B1,C1,D1],[A2,B2,C2,D2], [A3,B3,C3,D3], [A4,B4,C4,D4], Fun) ->
+ [Fun(A1,A2,A3,A4), Fun(B1,B2,B3,B4),
+ Fun(C1,C2,C3,C4), Fun(D1,D2,D3,D4)].
+
+apply6([A1,B1,C1,D1], [A2,B2,C2,D2],
+ [A3,B3,C3,D3], [A4,B4,C4,D4],
+ [A5,B5,C5,D5], [A6,B6,C6,D6],
+ Fun) ->
+ [Fun(A1,A2,A3,A4,A5,A6), Fun(B1,B2,B3,B4,B5,B6),
+ Fun(C1,C2,C3,C4,C5,C6), Fun(D1,D2,D3,D4,D5,D6)].
+
+apply7([A1,B1,C1,D1], [A2,B2,C2,D2],
+ [A3,B3,C3,D3], [A4,B4,C4,D4],
+ [A5,B5,C5,D5], [A6,B6,C6,D6], [A7,B7,C7,D7],
+ Fun) ->
+ [Fun(A1,A2,A3,A4,A5,A6,A7), Fun(B1,B2,B3,B4,B5,B6,B7),
+ Fun(C1,C2,C3,C4,C5,C6,C7), Fun(D1,D2,D3,D4,D5,D6,D7)].
+
+apply4c3([A1,B1,C1,D1], X,
+ [A3,B3,C3,D3], Y,
+ [A5,B5,C5,D5], Z,
+ [A7,B7,C7,D7],
+ Fun) ->
+ [Fun(A1,X,A3,Y,A5,Z,A7), Fun(B1,X,B3,Y,B5,Z,B7),
+ Fun(C1,X,C3,Y,C5,Z,C7), Fun(D1,X,D3,Y,D5,Z,D7)].
+
+
+get_node(Index, Tree) ->
+ Skip = (Index*?QNODE_SZ),
+ <<_:Skip/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,
+ _/binary>> = Tree,
+ #qnode{children=[C0,C1,C2,C3],
+ bboxes={[V00x,V10x,V20x,V30x],
+ [V00y,V10y,V20y,V30y],
+ [V00z,V10z,V20z,V30z],
+ [V01x,V11x,V21x,V31x],
+ [V01y,V11y,V21y,V31y],
+ [V01z,V11z,V21z,V31z]}}.
+
+get_qtriangles(Index, Binary) ->
+ Skip = (Index*?QTRI_SZ),
+ <<_:Skip/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,
+
+ E00x:?F32,E10x:?F32,E20x:?F32,E30x:?F32,
+ E00y:?F32,E10y:?F32,E20y:?F32,E30y:?F32,
+ E00z:?F32,E10z:?F32,E20z:?F32,E30z:?F32,
+
+ E01x:?F32,E11x:?F32,E21x:?F32,E31x:?F32,
+ E01y:?F32,E11y:?F32,E21y:?F32,E31y:?F32,
+ E01z:?F32,E11z:?F32,E21z:?F32,E31z:?F32,
+
+ Prim1:?I32, Prim2:?I32, Prim3:?I32, Prim4:?I32,
+ _/binary>> = Binary,
+ #qtri{v = {[V00x,V10x,V20x,V30x],
+ [V00y,V10y,V20y,V30y],
+ [V00z,V10z,V20z,V30z]},
+ e1 = {[E00x,E10x,E20x,E30x],
+ [E00y,E10y,E20y,E30y],
+ [E00z,E10z,E20z,E30z]},
+ e2 = {[E01x,E11x,E21x,E31x],
+ [E01y,E11y,E21y,E31y],
+ [E01z,E11z,E21z,E31z]},
+ f = [Prim1,Prim2, Prim3, Prim4]}.
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+print_tree(Nodes, PrimIndx, WorldBB) ->
+ print_tree(0, Nodes, PrimIndx, WorldBB, 1).
+
+print_tree(Index, Nodes, PrimIndx, Box, Level) ->
+ #qnode{children=Children, bboxes=BBs} = array:get(Index, Nodes),
+ [print_child(C, Nodes, PrimIndx, Box, 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, PBB, Level) ->
+ check_bb(PBB,BB),
+ 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, PBB, Level) ->
+ check_bb(PBB, BB),
+ io:format("~*.1c ~p ~s~n", [Level,$:, Node, f(BB)]),
+ print_tree(Node, Nodes, PrimIndx, BB, Level+2).
+
+check_bb(PBB,BB) -> check_bb(c, PBB,BB).
+check_bb(What, PBB, BB) -> %% Assert tree while printing it
+ case e3d_bv:union(PBB, BB) of
+ PBB -> ok;
+ _Other ->
+ io:format("Bad bounding box:~p ~n >>> ~s~n <<< ~s~n",
+ [What, f(PBB), f(BB)])
+ end.
+
+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("~.3g", [X]);
+f(X) when is_integer(X) ->
+ io_lib:format("~p", [X]);
+f(X) when is_list(X) ->
+ "[" ++ [f(E) ++ "," || E <- X] ++ "]".
+
+
+print_split(Start,Store,End, PrimIndx) ->
+ io:format("Left: [",[]),
+ [io:format("~p ", [array:get(Id,PrimIndx)])
+ || Id <- lists:seq(Start,Store-1)],
+ io:format("]~n Right: [",[]),
+ [io:format("~p ", [array:get(Id,PrimIndx)])
+ || Id <- lists:seq(Store,End-1)],
+ io:format("]~n",[]).
+
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
View
13 e3d/e3d_transform.erl
@@ -132,15 +132,16 @@ mul(#e3d_transf{mat=M1,inv=I1}, #e3d_transf{mat=M2,inv=I2}) ->
-spec lookat(e3d_point(), e3d_vector(), e3d_vector()) -> e3d_transform().
lookat(Pos, Look, Up) ->
Dir = e3d_vec:norm_sub(Look, Pos),
- Right = e3d_vec:norm(e3d_vec:cross(Dir, Up)),
+ Right = e3d_vec:norm(e3d_vec:cross(Dir, e3d_vec:norm(Up))),
NewUp = e3d_vec:norm(e3d_vec:cross(Right, Dir)),
- AsList = [tuple_to_list(Right), 0.0,
- tuple_to_list(NewUp), 0.0,
- tuple_to_list(Dir), 0.0,
- tuple_to_list(e3d_vec:neg(Pos)), 1.0],
+ AsList = [tuple_to_list(Right), 0.0,
+ tuple_to_list(NewUp), 0.0,
+ tuple_to_list(e3d_vec:neg(Dir)), 0.0,
+ 0.0, 0.0, 0.0, 1.0],
CamToWorld = list_to_tuple(lists:flatten(AsList)),
WorldToCam = e3d_mat:invert(CamToWorld),
- #e3d_transf{mat=WorldToCam,inv=CamToWorld}.
+ Translate = translate(identity(), e3d_vec:neg(Pos)),
+ mul(Translate, #e3d_transf{mat=WorldToCam,inv=CamToWorld}).
%%--------------------------------------------------------------------
%% @doc Generates a ortho transformation
View
63 e3d/test/qbvh.erl
@@ -9,7 +9,7 @@
-module(qbvh).
-compile(export_all).
--include_lib("e3d.hrl").
+-include_lib("../e3d.hrl").
-include_lib("cl/include/cl.hrl").
-define(F32, 32/float-native).
@@ -33,6 +33,12 @@
-define(RAYHIT_SZ, 16).
+-record(ray,
+ {o,d, % Origin, Direction vector
+ n, f}). % Near far (or MinT MaxT)
+
+-record(hit, {t, b1, b2, f = 16#ffffffff}).
+
go() -> start().
start() ->
@@ -42,15 +48,18 @@ start() ->
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)}
+ {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),
+ %% {Rays, Res} = init_cl(Qbvh),
+ %% check_rays(Rays, Res, array:from_list(Fs), ?VS),
+
+ erl_ray_trace([e3d_qbvh:ray({0.0, 0.5, 5.0},{0.0,0.0,-1.0})], Qbvh),
+ %%erl_ray_trace(Rays, Qbvh),
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) ->
+check_rays([#ray{d=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",
@@ -59,6 +68,20 @@ check_rays([Dir|Ds], <<T:?F32,B1:?F32,B2:?F32,Face:?I32, Rest/binary>>, Fs, Vs)
check_rays([],<<>>, _, _) ->
ok.
+erl_ray_trace([Ray|Rays], Qbvh) ->
+ case e3d_qbvh:ray_trace(Ray,Qbvh) of
+ #hit{f=16#ffffffff} ->
+ io:format("E miss~n",[]),
+ ok;
+ #hit{t=T, b1=B1, b2=B2, f=Face} ->
+ io:format("E ~s => ~s ~s ~s ~s~n",
+ [f(Ray#ray.d), f(T),f(B1),f(B2),f(Face)])
+ end,
+ erl_ray_trace(Rays, Qbvh);
+erl_ray_trace([], _) ->
+ 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))]),
@@ -100,7 +123,8 @@ init_cl(Qbvh) ->
end.
compile(CL, File) ->
- {ok, Bin} = file:read_file(File),
+ Dir = filename:dirname(code:which(?MODULE)),
+ {ok, Bin} = file:read_file(filename:join([Dir, File])),
{ok, Program} = clu:build_source(CL, Bin),
{ok, [Kernel]} = cl:create_kernels_in_program(Program),
Kernel.
@@ -128,7 +152,9 @@ create_rays() ->
MaxT = ?E3D_INFINITY,
Dirs = [calc_dir(Dir/255) || Dir <- lists:seq(0, 255)],
{<< << (create_ray(Origo, Dir, MinT, MaxT))/binary >>
- || Dir <- Dirs >>, Dirs, 256}.
+ || Dir <- Dirs >>,
+ [e3d_qbvh:ray(Origo, Dir, MinT, MaxT) || Dir <- 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}))).
@@ -137,3 +163,26 @@ 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>>.
+
+max_float() ->
+ max_float(3.000000e+38, 1.0e+37, 50000).
+
+max_float(X, Next, N) when N > 0, Next > 1.0 ->
+ Bin = <<X:?F32>>,
+ Op = try
+ <<Y:?F32>> = Bin,
+ erlang:display({Y, Next}),
+ inc
+ catch error:_ ->
+ erlang:display({fail,X,Next}),
+ dec
+ end,
+ if
+ Op == inc ->
+ max_float(X+Next, Next*1.2, N);
+ Op == dec ->
+ max_float(X-Next, Next / 2, N-1)
+ end;
+max_float(X,_,_) ->
+ X.
+
Please sign in to comment.
Something went wrong with that request. Please try again.