function I = sloray
% SLORAY Traçado de raio básico para fins educativos.
% output:
%   I: imagem RGB (matriz double MxNx3 com valores entre 0-1)
% Referências:
% [1] K. Suffern. Ray Tracing from the Ground Up. 2007.
% [2] A. Watt, M. Watt. Advanced Animation and Rendering Techniques. 1992.
% [3] P. Heckbert. Writing a Ray Tracer. In: A. Glassner, An Introduction to Ray Tracing. 1989

% EXEMPLO DE USO:
% tic, I=sloray; toc, figure, imshow(I)

% SloRay - Slow Raytracer
% Autores:
% Augusto Valente (avalente arr_oba dca ponto fee ponto unicamp ponto br)
% Wallace Loos (wloos arr_oba dca ponto fee ponto unicamp ponto br)
% Junho/2011

% tamanho da imagem em pixels
% - largura
res_w = 128;
% - altura
res_h = res_w*3/4; % (usando aspect ratio 1.333)
% tamanho do plano de imagem centrado na origem (coords. do "mundo")
% - largura
win_w = 4;
% - altura
win_h = win_w*3/4;

% posicao do observador
ray_o = [0 0 2];
% "epsilon factor": [1], 16.4
kEps = 1e-5;
% iluminacao ambiente ([2], cap. 2)
Ia_ka = [.2 .2 .2];
% máximo nível de recursão
MAX_DEPTH = 5;
% mínimo peso de um raio para levá-lo em conta na recursão
MIN_WEIGHT = 1e-3;
% cor de fundo
BACKGROUND_COLOR = [0 0 0];
% estatística
stat_max_depth = 0;
% posição no eixo z do plano de imagem
Z_IMAGE_PLANE = 0;

% posicoes dos pixels no plano de imagem
pixel_w = win_w/res_w;
pixel_h = win_h/res_h;
c_x = -win_w/2+pixel_w/2:pixel_w:win_w/2-pixel_w/2;
c_y = win_h/2-pixel_h/2:-pixel_h:-win_h/2+pixel_h/2;
% c_x = linspace(-win_w/2,win_w/2,res_w);
% c_y = linspace(win_h/2,-win_h/2,res_h);

% imagem
I = zeros(res_h, res_w, 3);

% estrutura "objs" com objetos da cena
% campos: 
% - "parms": cell-array com 2 parametros geometricos
%   --> esfera: {centro raio}
%   --> plano:  {ponto normal}
% - "cor": função que recebe uma posição e retorna uma cor RGB (usado como coeficiente kd do modelo de Phong)
% - "krl": coeficiente de especularidade do modelo de Phong (local)
% - "ns":  coeficiente de concentração da especularidade do modelo de Phong (local)
% - "krg": coeficiente de especularidade global
% - "intersec": função que recebe parâmetros do raio e do objeto e testa a intersecção

% 3 objetos: 
% esfera especular azul (R = 1, centro = (-2, 0, 4))
% plano lambertiano com textura xadrez (normal = (.5, 1, 0), ponto = (0, -3, -4))
% esfera lambertiana verde (R = 0.5, centro = (1, 0, -4))
objs = struct( ...
    'parms', {{[-2 0 -4] 1}, {[0 -3 -4] [.5 1 0]/norm([.5 1 0])}, {[1 0 -4] .5}}, ...
    'cor', {@cor1, @checker, @cor2}, ...
    'krl', {.9, 0, 0}, ...
    'ns', {40, 0, 0}, ...
    'krg', {.3, 0, 0}, ...
    'intersec', {@intersec_sphere, @intersec_plane, @intersec_sphere}, ...
    'normal', {@normal_sphere, @normal_plane, @normal_sphere});

% estrutura com fontes de luz pontuais
% campos:
% - "pos": posição
% - "I": intensidade RGB
lights = struct('pos', {[10 10 -4]},'I',{[245,255,250]/255});
% lights = struct('pos', {[10 10 -4] [0 20 -8]},'I',{[245,255,250]/255,[245,255,250]/255});


h_wait = waitbar(0,'aguarde...');
% traça 1 raio por pixel
for ii=1:res_w
for jj=1:res_h
    % direção do raio
    ray_d = [c_x(ii) c_y(jj) Z_IMAGE_PLANE]-ray_o;
    ray_d = ray_d / norm(ray_d);
    % começa recursão e seta pixel da imagem
    I(jj,ii,:) = trace_ray(ray_d, ray_o, 1, 1);
end
waitbar(ii/res_w,h_wait)
end
close(h_wait)
disp(['Nível de recursão máximo: ' num2str(stat_max_depth) '/' num2str(MAX_DEPTH)])

function cor = trace_ray(ray_d, ray_o, depth, weight)
% TRACE_RAY Traça um raio e determina sua contribuição em cor.
% input:
%   ray_d: direção do raio (vetor 1x3)
%   ray_o: origem do raio (vetor 1x3)
%   depth: nível de recursão
%   weight: fator de atenuação da energia após reflexões
% output:
%   cor: cor RGB (vetor 1x3)

if depth > MAX_DEPTH || weight < MIN_WEIGHT, cor = BACKGROUND_COLOR; return; end
if depth > stat_max_depth, stat_max_depth = depth; end

% pega hit mais perto e faz shading
[obj_hit, t] = ray_hit(ray_d, ray_o);
if obj_hit
    cor = shade(ray_d, ray_o, obj_hit, t, depth, weight);
else
    cor = BACKGROUND_COLOR;
end
end

function [obj_hit, t] = ray_hit(ray_d, ray_o)
% RAY_HIT Testa interseções para um dado raio
%         com todos os objetos da struct global 'objs'.
% input:
%   ray_d: direção do raio (vetor 1x3)
%   ray_o: origem do raio (vetor 1x3)
% output:
%   obj_hit: índice do objeto interceptado (0 c.c.)
%   t: satisfaz a equação P_intersec = ray_o + t * ray_d

obj_hit = 0;
t = inf;
for k=1:length(objs)
    [hit, t1] = objs(k).intersec(ray_d, ray_o, objs(k).parms);
    if hit && t1 < t
        t = t1;
        obj_hit = k;
    end
end
end

function [hit, t] = intersec_sphere(ray_d, ray_o, obj_parms)
% INTERSEC_SPHERE Interseção raio-esfera. ref: [1], p. 58
% input:
%   ray_d: direção do raio (vetor 1x3)
%   ray_o: origem do raio (vetor 1x3)
%   obj_parms: cell-array 1x2 contendo centro (vetor 1x3) e raio (escalar)
% output:
%   hit: flag de interseção
%   t: satisfaz a equação P_intersec = ray_o + t * ray_d

obj_c = obj_parms{1};
obj_r = obj_parms{2};

a = dot(ray_d,ray_d);
b = 2*dot(ray_o - obj_c, ray_d);
c = dot( ray_o - obj_c, ray_o - obj_c ) - obj_r^2;
delta = b^2-4*a*c;

t = inf;
hit = false;
if delta >= 0
    t = (-b-sqrt(delta)) / (2*a);
    if t > kEps
        hit = true;
        return
    end
    t = (-b+sqrt(delta)) / (2*a);
    if t > kEps
        hit = true;
        return
    end
end
end

function [hit, t] = intersec_plane(ray_d, ray_o, obj_parms)
% INTERSEC_PLANE Interseção raio-plano. ref: [1], p. 56
% input:
%   ray_d: direção do raio (vetor 1x3)
%   ray_o: origem do raio (vetor 1x3)
%   obj_parms: cell-array 1x2 contendo um ponto (vetor 1x3) e a normal
%       (vetor 1x3) do plano
% output:
%   hit: flag de interseção
%   t: satisfaz a equação P_intersec = ray_o + t * ray_d

a = obj_parms{1};
n = obj_parms{2};
t = dot(a - ray_o, n) / dot(ray_d, n);
hit = false;
if t > kEps
    hit = true;
end
end

function cor = cor1(~)
% COR1 Gera verde RGB.
cor = [0 1 0];
end

function cor = cor2(~)
% COR1 Gera azul RGB.
cor = [0 0 1];
end

function cor = checker(P)
% CHECKER Gera textura quadriculada usando floor(). ref: [1], p. 674
% input:
%   P: ponto no espaço (vetor 1x3)
% output:
%   cor: cor RGB (vetor 1x3)

P = P - 0.000187453738; % evita que o plano caia em região flat
size = .2; % lado dos quadrados
if mod(sum(floor(P/size)),2) == 0, cor = [139,90,0]/255;
else cor = 	[238,154,0]/255;
end
end

function cor = shade(ray_d, ray_o, obj_hit, t, depth, weight)
% SHADE Calcula cor em uma interseção com objeto.
% input:
%   ray_d: direção do raio (vetor 1x3)
%   ray_o: origem do raio (vetor 1x3)
%   obj_hit: índice do objeto em questão
%   t: satisfaz a equação P_intersec = ray_o + t * ray_d
%   depth: nível de recursão
%   weight: fator de atenuação da energia após reflexões
% output:
%   cor: cor RGB (vetor 1x3)

obj = objs(obj_hit);
P = ray_o + t * ray_d; % ponto de interseção
n = obj.normal(obj, P); % normal à superfície
cor = Ia_ka.*obj.cor(P);
% determina contribuição de cada fonte de luz
for l=1:length(lights)
    L = lights(l).pos - P; % shadow feeler
    dist = norm(L); % distância à fonte de luz
    dot_L_n = dot(L/dist,n);
    % verifica se há oclusão
    [obj_hit, t] = ray_hit(L/dist, P);
    if (obj_hit && t < dist) % ver [1], fig. 16.8
        continue;
    end
    if dot_L_n > 0
        % modelo de Phong
        cor = cor + lights(l).I .* obj.cor(P) * dot_L_n; % [2], eq. 2.5
    end
    % direção especular utilizada no modelo de Phong
    R_light = 2*dot_L_n*n - L/dist;
    dot_R_V = dot(R_light,-ray_d);
    if dot_R_V > 0
        % modelo de Phong
        cor = cor + lights(l).I .* obj.krl * dot_R_V^obj.ns; % [2], eq. 2.5
    end
end
% direção do raio refletido
R = ray_d - 2*dot(ray_d,n)*n;
% contribuição do raio refletido (recursão)
cor = cor + weight*obj.krg*trace_ray(R, P, depth+1, weight*obj.krg);
end

function n = normal_sphere(obj, P)
% NORMAL_SPHERE Normal à superfície da esfera no ponto.
% input:
%   obj: item da struct global 'objs'
%   P: ponto na superfície da esfera (vetor 1x3)
% output:
%   n: normal
% parâmetros da estrutura:
%   obj.parms{1}: centro da esfera (vetor 1x3)
%   obj.parms{2}: raio da esfera

% n = (P - C) / R (vetor unitário do centro C para o ponto P)
n = (P - obj.parms{1}) / obj.parms{2};
end

function n = normal_plane(obj, ~)
% NORMAL_PLANE Normal ao plano.
% input:
%   obj: item da struct global 'objs'
% output:
%   n: normal
% parâmetros da estrutura:
%   obj.parms{2}: normal unitária ao plano (vetor 1x3)
n = obj.parms{2};
end

end