レイトレースベンチ完成

シャドウレイの判定を作っているところで、いろいろ間違っている部分が発覚して直すのに手間取ってしまった。
後ろ向きにレイートレースしていたりとか、
画素の順番がおかしかったりとかいろいろな罠にはまった。
リスト内包表記で複数のリストを組み合わせるときの順番は制御できないのだろうか。
終盤は座標を確認しながら丹念にプリントデバッグすることになったけど、実行時エラーで行番号が判らないのは厄介だった。
エラー処理も勉強しておく必要があるなぁ。

次の課題は、これを並列化して更には分散化すること。
erlang的にはお楽しみタイムですｗ

-module(renderer).
-export([compatible/3]).

% export
%------------------------------------------------------------%
% レイトレースベンチ互換の動作
compatible(F, Size, Level)->
  Scene=create_scene(Level),
  %io:write(Scene), io:format("~n"),
  %io:write(create_screen(Size, Size)), io:format("~n"),
  dib32:write_file(F, Size, Size, 
      lists:map(fun(PIXEL) -> sampling(PIXEL, Scene, 4, ray_generator(Size, Size, Size)) end, 
        create_screen(Size, Size))).

% サンプリング
%------------------------------------------------------------%
% 幅と高さを指定してピクセル座標の配列を取得する
create_screen(W, H)->
  lists:sort(fun({X1, Y1}, {X2, Y2}) -> 
    if Y1<Y2 -> true;
       Y1>Y2 -> false;
       true -> X1<X2 
     end
       end,
  [{X, Y}||X<-generate_list(W), Y<-generate_list(H)]).
  

% 指定ピクセルをサンプリング
% ピクセルの中央をサンプリング
sampling({PIXEL_X, PIXEL_Y}, Scene, 1, GET_RAY)->
  ray_trace(GET_RAY(PIXEL_X+0.5, PIXEL_Y+0.5), Scene);
% オーバーサンプリング(N x N倍サンプリングして平均する)
sampling({PIXEL_X, PIXEL_Y}, Scene, N, GET_RAY)->
  average_color(
    lists:map(fun({X, Y})->ray_trace(GET_RAY(X, Y), Scene) end, 
      over_sampling({PIXEL_X, PIXEL_Y}, N))
  ).

% オーバーサンプリング位置を生成
over_sampling({PIXEL_X, PIXEL_Y}, N)->
  D=1.0/N,
  [{X, Y}|| X<-append_times(N, D, [PIXEL_X]), Y<-append_times(N, D, [PIXEL_Y])].

% スクリーン位置XYに対応するレイを作成する関数を返す
ray_generator(W, H, Z)->
  fun (X, Y)->
    [
    {orig, [0.0, 0.0, -4.0]}, % 視点
    {dir, vec:normalize([X-W/2, Y-H/2, Z])} % 方向
    ]
  end.

% 指定されたピクセルにレイを飛ばしRGB値を得る
ray_trace(Ray, Scene)->
  Intersect=intersect(Ray, get_value(Scene, objects)),
  %io:write(PIXEL), io:write(Intersect), io:format("~n"),
  case Intersect of
    % 交差せず
    {false} -> 
      [0, 0, 0, 0];
    % 交差した
    _ -> 
      shading(Intersect, Scene)
  end.

% シーン
%------------------------------------------------------------%
% シーンを作成する
create_scene(N)->
  [
  {lights, [ {directional_light, [{dir, vec:normalize([-1.0, -3.0, 2.0])}]} ]} , 
  {objects, create_spheres(N, [0, -1, 1], 1, [])}
  ].

% レベルが1になるまで再帰的に球を作る
create_spheres(N, [X, Y, Z], R, L) when N>1 ->
  % 中央
  L2=[create_sphere([X, Y, Z], R) | L],
  % ずれる距離
  F=3*R/math:sqrt(12),
  L3=create_spheres(N-1, [X+F, Y+F, Z+F], R/2, L2),
  L4=create_spheres(N-1, [X-F, Y+F, Z+F], R/2, L3),
  L5=create_spheres(N-1, [X-F, Y+F, Z-F], R/2, L4),
  create_spheres(N-1, [X+F, Y+F, Z-F], R/2, L5);
% 再帰終了
create_spheres(1, P, R, L)->
  [create_sphere(P, R)|L].

% 球を作成
create_sphere(P, R)->
  {sphere, [{center, P}, {radius, R}]}.

% 交差判定
%------------------------------------------------------------%
intersect(Ray, [H|T], Result)->
  %io:format("Ray:"), io:write(Ray), io:format("~n"),
  [intersect(Ray, H)|intersect(Ray, T, Result)];
intersect(_, [], Result)->
  Result.

% 球の交差判定
intersect(Ray, {sphere, Sphere})->
   %io:write(Ray), io:format("Sphere~n"),
   % 視点から球の中心へのベクトル
   Vec=vec:sub(get_value(Sphere, center), get_value(Ray, orig)),
   B=vec:dot(Vec, get_value(Ray, dir)),
   R=get_value(Sphere, radius),
   case B*B-vec:sqnorm(Vec)+R*R of
     SqDet when SqDet<0 -> 
       {false}; % 交差しない
     SqDet ->
       % 交差したので交点の距離と法線を求める
       Det=math:sqrt(SqDet),
       T2=B+Det,
       if 
         T2<0 -> {false};
         true ->
           T1=B-Det,
           if 
             T1>0 ->
               % T1採用
               T=T1;
             true ->
               % T2採用
               T=T2
           end,
           % 交点を求める
           P=vec:add(get_value(Ray, orig), vec:times(T, get_value(Ray, dir))),
           % 法線は交点-球の中心
           N=vec:normalize(vec:sub(P, get_value(Sphere, center))),
           {hit, [
            {distance, T}, 
            {normal, N},
            {pos, vec:add(P, vec:times(0.001, N))} % 法線方向に浮かせた交点
            %{pos, P}
          ]}
      end
   end;
% リストの交差判定
intersect(Ray, L)->
  case [X||{hit, X}<-intersect(Ray, L, [])] of
    [] -> {false};
    Intersection->
      lists:nth(1, lists:sort(fun(A, B)-> 
        get_value(A, distance)<get_value(B, distance) end,
        Intersection))
  end.

% 輝度計算
%------------------------------------------------------------%
shading(Intersection, Scene)->
  shading(Intersection, Scene, get_value(Scene, lights), [0, 0, 0, 0]).
% リスト積算
shading(Intersection, Scene, [Light|T], Diffuse)->
  shading(Intersection, Scene, T, vec:add(Diffuse, calc_shading(Intersection, Scene, Light)));
% リスト終端
shading(_, _, [], Diffuse)->
  Diffuse.

calc_shading(Intersection, Scene, {directional_light, Light}) ->
  Cos=trunc(-255 * vec:dot(get_value(Light, dir), get_value(Intersection, normal))),
  if Cos<0 -> [0, 0, 0, 0] ;
     true -> 
       % 影判定
       ShadowRay=[
           {orig, get_value(Intersection, pos)}, % 交点から
           {dir, vec:negative(get_value(Light, dir))} % 光源に向かう
           ],
       %io:write(Intersection), io:format("~n"),
       %io:write(ShadowRay), io:format("~n"),
       Shadow=intersect(ShadowRay, get_value(Scene, objects)),
       %io:write(Shadow), io:format("~n"),
       case Shadow of
         % 光源に到達。影にならない
         {false} -> 
           [Cos, Cos, Cos, 255];
         % 遮られた。影
         _ -> 
           [0, 0, 0, 0]
       end
  end.

% その他
%------------------------------------------------------------%
% {atom, Param}のリストからatomの一致するものをゲットする
get_value([], _)->
  false;
get_value([H|T], Key)->
  {K, V}=H,
  if 
    K==Key -> V;
    true -> get_value(T, Key)
  end.

% 0からN-1のリストを作る
generate_list(N)->
  generate_list(N-1, []).
generate_list(0, L)->
  [0|L];
generate_list(N, L)->
  generate_list(N-1, [N|L]).

% 指定回数Dを加算した配列をつくる
append_times(0, _, L)->
  L;
append_times(N, D, [H|T])->
  append_times(N-1, D, [H+D|[H|T]]).

% RGBAを積算して平均を得る
average_color(L)->
  average_color(L, [0, 0, 0, 0], 0).
average_color([[R, G, B, A]|T], [SUMR, SUMG, SUMB, SUMA], N)-> 
  average_color(T, [SUMR+R, SUMG+G, SUMB+B, SUMA+A], N+1);
average_color([], [SUMR, SUMG, SUMB, SUMA], N)->
  [SUMR div N, SUMG div N, SUMB div N, SUMA div N].

警告が出ないように直していたらRGBAの平均のとこにタイポがあったので修正ｗ

227> renderer:compatible("tmp.bmp", 200, 3).
ok

順調に数分かかるようになりました。