from Vectorized raytracer for sphere objects by Brian Jensen
Vectorized raytracer for sphere objects

ray_tracer_final.m
%Raytracer for one sphere object and local illumination
%
%Made by Brian Jensen
%Aalborg University (Denmark) - Computer Vision and Graphics
%E-mail: smoller_@hotmail.com

clear;close all
tic

image_pixels_x = 2^(10)      %pixels in horisontal direction          (value: int)
image_pixels_y = 2^(10)      %pixels in vertical direction            (value: int)
%intersections is a memory and speed optimization(speedup depend on the imagesize)
intersections = 32;          %number of small subrender areas         (value: 2^n)


%objects 1 propertise
sphere_center = [-2 -2 20 1]';  %sphere center, radius, number

%object 1 materiale propertise (Left handed coordinate system)
material_diffuse = [1 0.0 0]';      %diffuse reflection
material_specular = [1 1 1]';       %specular reflection
material_ambient = [1 1 1]';        %ambient reflection
reflection_coef = 50;               %specular reflection coefficient

%lights 1 propertise
light_diffuse = [1 1 1]';           %diffuse color
light_specular = [1 1 1]';          %specular color
light_ambient = [0.3 0.2 0.2]';     %ambient color
light_pos = sphere_center(1:3,1) + [-2 10 -10]';  %light position relativ to object 1

%view plane and eye
center_point = [0 0 0]';            %viewplane centerpoint
window_ratio = normc([image_pixels_x/image_pixels_y 1]');   %viewplane ratio
window_scale = [1/image_pixels_x 1/image_pixels_y];         %viewplane scale
window_dist = 1;                    %viewplane distance raletiv to z axis for center point

area_render_x = image_pixels_x / intersections;         %small viewplane area x
area_render_y = image_pixels_y / intersections;         %small viewplane area y

render_image = zeros(image_pixels_y,image_pixels_x,3);  %init image array

%iterations for all small viewplanes
for area_x_iteration = 1 : intersections
    %iterations for all small viewplanes
    for area_y_iteration = 1 : intersections

%viewplane stuff - start
pixel_coordinate_x = ones(1,area_render_y)'*[1+(area_x_iteration-1)*area_render_x:area_x_iteration*area_render_x];
pixel_coordinate_y = [1+(area_y_iteration-1)*area_render_y : area_y_iteration*area_render_y]'*ones(1,area_render_x);

window_pixel(:,:,1) = [window_ratio(1)*(pixel_coordinate_x - image_pixels_x/2)*window_scale(1)];
window_pixel(:,:,2) = [window_ratio(2)*(pixel_coordinate_y - image_pixels_y/2)*window_scale(2)];
window_pixel(:,:,3) = [window_dist*ones(area_render_y,area_render_x)];
%viewplane stuff - end

%object detection / distance - start
%should be possible to optimize - start
for i = 1:3
    delta(:,:,i) = window_pixel(:,:,i) - center_point(i); %delta
end
A = sum(delta.^2,3);

for i = 1:3
    B(:,:,i) = ( delta(:,:,i)*(center_point(i)-sphere_center(i,1)) );
end
B = 2*sum(B,3);
%should be possible to optimize - end

C = (sum( (center_point-sphere_center(1:3,1)).^2 ) - sphere_center(4,1)^2);
D = B.^2 -(4*A.*C);

dist = (-B+D.^0.5) ./ (2.*A+eps);
isreal_test = (~(imag(dist) ~= 0));
dist = dist.*isreal_test;
%object detection / distance - end

%local illumination - start
reshaped_delta = reshape(delta,[area_render_x*area_render_y 3]);
reshape_array = ones(1,(area_render_x*area_render_y));

ray_hit_coordinate = reshaped_delta.*(  reshape(dist ,[area_render_x*area_render_y 1])*ones(1,3) );
v_norm = normr(ray_hit_coordinate - reshaped_delta);
n_norm = normr(ray_hit_coordinate - reshape_array'*(sphere_center(1:3,1) - center_point)' );
l_norm = normr(ray_hit_coordinate - reshape_array'*(light_pos - center_point)' );
h_norm = normr(v_norm + l_norm);
     
reshaped_Ia = reshape_array'*(material_ambient .* light_ambient)';
reshaped_Id = reshape_array'*(material_diffuse .* light_diffuse)' .*  ((max([dot(n_norm,l_norm,2) zeros(area_render_x*area_render_y,1)],[],2))*ones(1,3));
reshaped_Is = reshape_array'*(material_specular .* light_specular)' .*  ((max([dot(n_norm,h_norm,2) zeros(area_render_x*area_render_y,1)],[],2).^reflection_coef)*ones(1,3));
reshaped_Ir = reshaped_Ia + reshaped_Id + reshaped_Is;
    
isreal_tests(:,:,1) = isreal_test;
isreal_tests(:,:,2) = isreal_test;
isreal_tests(:,:,3) = isreal_test;
Ir = reshape(reshaped_Ir,[area_render_y area_render_x 3]).*isreal_tests;
%local illumination - end

%add small viewplane to image array
render_image(1+(area_y_iteration-1)*area_render_y : area_y_iteration*area_render_y,1+(area_x_iteration-1)*area_render_x:area_x_iteration*area_render_x,:) = Ir;

    end
end
toc

%imagesc(isreal_test)
image(render_image./max(max(max(render_image))))
axis equal
%imwrite(render_image,'c:\test_render.bmp') 

Contact us at files@mathworks.com