% Further testing for building a radiosity lighting method for virtual
% spaces. Here we pixelize all the walls and draw them with a smooth
% gradient of gray values.
%
% Samuli Siltanen December 2012

% Choose mosaic resolution
n = 10;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Construct centerpoints of the mosaic tiles.
% The letter d denotes the length of the side of a pixel.
d        = 2/n;
tmp      = -1-d/2 + [1:n]*d;

% Initialize centerpoint coordinate matrices
Xmat = zeros(n^2,5);
Ymat = zeros(n^2,5);
Zmat = zeros(n^2,5);

% The back wall
[X,Z]     = meshgrid(tmp);
Xmat(:,1) = X(:);
Zmat(:,1) = Z(:);
Ymat(:,1) = ones(n^2,1);

% Roof
[X,Y] = meshgrid(tmp);
Xmat(:,2) = X(:);
Ymat(:,2) = Y(:);
Zmat(:,2) = ones(n^2,1);

% Floor
Xmat(:,3) = X(:);
Ymat(:,3) = Y(:);
Zmat(:,3) = -ones(n^2,1);

% Right-hand-side wall
[Y,Z] = meshgrid(tmp);
Ymat(:,4) = Y(:);
Zmat(:,4) = Z(:);
Xmat(:,4) = ones(n^2,1);

% Left-hand-side wall
Ymat(:,5) = Y(:);
Zmat(:,5) = Z(:);
Xmat(:,5) = -ones(n^2,1);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Construct the color matrix using the radiosity lighting model.

% Scaling parameter for the inverse square attenuation law. 
% Taking sc_par = 1; leads to unnatural lighting, taking sc_par > 2 has
% better scaling between the size of the cubic room and the attenuation
% law.
sc_par = 3;

% Form the geometrical view factor matrix F. Note that F is symmetric. 
% See http://en.wikipedia.org/wiki/View_factor for details of computation.
F = zeros(5*n^2);
% From the roof (jjj) to the back wall (iii)
for iii = 1:n^2
    for jjj = 1:n^2
        % Centerpoint of the current pixel in the back wall
        piii = [Xmat(iii,1);Ymat(iii,1);Zmat(iii,1)];
        % Centerpoint of the current pixel in the roof
        pjjj = [Xmat(jjj,2);Ymat(jjj,2);Zmat(jjj,2)];
        % Distance between the points
        difvec = piii-pjjj;
        r      = norm(difvec);
        % View angles
        tmp2   = difvec/r;
        cosjjj = abs(tmp2(3));
        cosiii = abs(tmp2(2));
        % Calculate element of F
        F(iii,n^2+jjj) = cosiii*cosjjj/(pi*(sc_par*r)^2);
    end
end
% From the floor (jjj) to the back wall (iii)
for iii = 1:n^2
    for jjj = 1:n^2
        % Centerpoint of the current pixel in the back wall
        piii = [Xmat(iii,1);Ymat(iii,1);Zmat(iii,1)];
        % Centerpoint of the current pixel in the roof
        pjjj = [Xmat(jjj,3);Ymat(jjj,3);Zmat(jjj,3)];
        % Distance between the points
        difvec = piii-pjjj;
        r      = norm(difvec);
        % View angles
        tmp2   = difvec/r;
        cosjjj = abs(tmp2(3));
        cosiii = abs(tmp2(2));
        % Calculate element of F
        F(iii,2*n^2+jjj) = cosiii*cosjjj/(pi*(sc_par*r)^2);
    end
end
% From the right-hand-side wall (jjj) to the back wall (iii)
for iii = 1:n^2
    for jjj = 1:n^2
        % Centerpoint of the current pixel in the back wall
        piii = [Xmat(iii,1);Ymat(iii,1);Zmat(iii,1)];
        % Centerpoint of the current pixel in the roof
        pjjj = [Xmat(jjj,4);Ymat(jjj,4);Zmat(jjj,4)];
        % Distance between the points
        difvec = piii-pjjj;
        r      = norm(difvec);
        % View angles
        tmp2   = difvec/r;
        cosjjj = abs(tmp2(1));
        cosiii = abs(tmp2(2));
        % Calculate element of F
        F(iii,3*n^2+jjj) = cosiii*cosjjj/(pi*(sc_par*r)^2);
    end
end
% From the left-hand-side wall (jjj) to the back wall (iii)
for iii = 1:n^2
    for jjj = 1:n^2
        % Centerpoint of the current pixel in the back wall
        piii = [Xmat(iii,1);Ymat(iii,1);Zmat(iii,1)];
        % Centerpoint of the current pixel in the roof
        pjjj = [Xmat(jjj,5);Ymat(jjj,5);Zmat(jjj,5)];
        % Distance between the points
        difvec = piii-pjjj;
        r      = norm(difvec);
        % View angles
        tmp2   = difvec/r;
        cosjjj = abs(tmp2(1));
        cosiii = abs(tmp2(2));
        % Calculate element of F
        F(iii,4*n^2+jjj) = cosiii*cosjjj/(pi*(sc_par*r)^2);
    end
end
% From the floor (jjj) to the roof (iii)
for iii = 1:n^2
    for jjj = 1:n^2
        % Centerpoint of the current pixel in the back wall
        piii = [Xmat(iii,2);Ymat(iii,2);Zmat(iii,2)];
        % Centerpoint of the current pixel in the roof
        pjjj = [Xmat(jjj,3);Ymat(jjj,3);Zmat(jjj,3)];
        % Distance between the points
        difvec = piii-pjjj;
        r      = norm(difvec);
        % View angles
        tmp2   = difvec/r;
        cosjjj = abs(tmp2(3));
        cosiii = abs(tmp2(3));
        % Calculate element of F
        F(n^2+iii,2*n^2+jjj) = cosiii*cosjjj/(pi*(sc_par*r)^2);
    end
end
% From the right-hand wall (jjj) to the roof (iii)
for iii = 1:n^2
    for jjj = 1:n^2
        % Centerpoint of the current pixel in the back wall
        piii = [Xmat(iii,2);Ymat(iii,2);Zmat(iii,2)];
        % Centerpoint of the current pixel in the roof
        pjjj = [Xmat(jjj,4);Ymat(jjj,4);Zmat(jjj,4)];
        % Distance between the points
        difvec = piii-pjjj;
        r      = norm(difvec);
        % View angles
        tmp2   = difvec/r;
        cosjjj = abs(tmp2(1));
        cosiii = abs(tmp2(3));
        % Calculate element of F
        F(n^2+iii,3*n^2+jjj) = cosiii*cosjjj/(pi*(sc_par*r)^2);
    end
end
% From the left-hand-side wall (jjj) to the roof (iii)
for iii = 1:n^2
    for jjj = 1:n^2
        % Centerpoint of the current pixel in the back wall
        piii = [Xmat(iii,2);Ymat(iii,2);Zmat(iii,2)];
        % Centerpoint of the current pixel in the roof
        pjjj = [Xmat(jjj,5);Ymat(jjj,5);Zmat(jjj,5)];
        % Distance between the points
        difvec = piii-pjjj;
        r      = norm(difvec);
        % View angles
        tmp2   = difvec/r;
        cosjjj = abs(tmp2(1));
        cosiii = abs(tmp2(3));
        % Calculate element of F
        F(n^2+iii,4*n^2+jjj) = cosiii*cosjjj/(pi*(sc_par*r)^2);
    end
end
% From the right-hand-side wall (jjj) to the floor (iii)
for iii = 1:n^2
    for jjj = 1:n^2
        % Centerpoint of the current pixel in the back wall
        piii = [Xmat(iii,3);Ymat(iii,3);Zmat(iii,3)];
        % Centerpoint of the current pixel in the roof
        pjjj = [Xmat(jjj,4);Ymat(jjj,4);Zmat(jjj,4)];
        % Distance between the points
        difvec = piii-pjjj;
        r      = norm(difvec);
        % View angles
        tmp2   = difvec/r;
        cosjjj = abs(tmp2(1));
        cosiii = abs(tmp2(3));
        % Calculate element of F
        F(2*n^2+iii,3*n^2+jjj) = cosiii*cosjjj/(pi*(sc_par*r)^2);
    end
end
% From the left-hand-side wall (jjj) to the floor (iii)
for iii = 1:n^2
    for jjj = 1:n^2
        % Centerpoint of the current pixel in the back wall
        piii = [Xmat(iii,3);Ymat(iii,3);Zmat(iii,3)];
        % Centerpoint of the current pixel in the roof
        pjjj = [Xmat(jjj,5);Ymat(jjj,5);Zmat(jjj,5)];
        % Distance between the points
        difvec = piii-pjjj;
        r      = norm(difvec);
        % View angles
        tmp2   = difvec/r;
        cosjjj = abs(tmp2(1));
        cosiii = abs(tmp2(3));
        % Calculate element of F
        F(2*n^2+iii,4*n^2+jjj) = cosiii*cosjjj/(pi*(sc_par*r)^2);
    end
end
% From the left-hand-side wall (jjj) to the right-hand-side wall (iii)
for iii = 1:n^2
    for jjj = 1:n^2
        % Centerpoint of the current pixel in the back wall
        piii = [Xmat(iii,4);Ymat(iii,4);Zmat(iii,4)];
        % Centerpoint of the current pixel in the roof
        pjjj = [Xmat(jjj,5);Ymat(jjj,5);Zmat(jjj,5)];
        % Distance between the points
        difvec = piii-pjjj;
        r      = norm(difvec);
        % View angles
        tmp2   = difvec/r;
        cosjjj = abs(tmp2(1));
        cosiii = abs(tmp2(1));
        % Calculate element of F
        F(3*n^2+iii,4*n^2+jjj) = cosiii*cosjjj/(pi*(sc_par*r)^2);
    end
end

% Add the contribution of the area of each pixel
F = (d^2)*F;

% Use symmetry to finish the construction of F
F = F+F.';

% Construct the right hand side Evec of the radiosity equation. Evec
% describes the contribution of emitted light in the scene. For example,
% each pixel belonging to a lamp in the virtual space causes a positive 
% element in Evec.
Evec   = zeros(5*n^2,1);
indvec = repmat(logical(0),size(Evec));
indvec(n^2+[1:n^2]) = sqrt((Xmat(:,2)-.3).^2+Ymat(:,2).^2)<.3;
Evec(indvec) = 1;
% Evec(n^2+round(n^2/2)-2) = 1;
% Evec(3*n^2+round(n^2/2)-2) = 1;

% Solve for color vector
rho = 1;
colorvec = (eye(5*n^2)-rho*F)\Evec;

% Normalize the values of the color vector between 0 and 1.
% The cut_param below is for preventing the lamp color being too far away
% from the gray colors of everything else. Setting cut_param = 1; has no
% effect, and setting 0 < cut_param < 1 can improve the plot. Experiment 
% with it to find the best value.
colorvec  = max(0,colorvec);
cut_param = .5;
colorvec  = colorvec/(cut_param*max(colorvec));
colorvec  = min(1,colorvec);

% Gamma correction for optimal gray levels. Choosing gammacorr = 1; has no
% effect, and taking 0 < gammacorr = < 1 small will lighten up the dark
% shades.
gammacorr = .25;
colorvec  = colorvec.^(gammacorr);

% Construct color matrix by repetition
colormat = repmat(colorvec(:),1,3);





%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Draw the illuminated virtual space

% Create plot window
figure(1)
clf

% Draw all the walls consisting of n x n little squares (pixels). 
% Pick the gray value of each square from the illumination vector
% calculated by the radiosity method above

% The back wall
colorind = 1;
for iii = 1:(n^2)
    patch(...
        [Xmat(iii,1)+d/2,Xmat(iii,1)+d/2,Xmat(iii,1)-d/2,Xmat(iii,1)-d/2],...
        [Ymat(iii,1),Ymat(iii,1),Ymat(iii,1),Ymat(iii,1)],...
        [Zmat(iii,1)-d/2,Zmat(iii,1)+d/2,Zmat(iii,1)+d/2,Zmat(iii,1)-d/2],...
        colormat(colorind,:))
    colorind = colorind+1;
end

% Roof
for iii = 1:(n^2)
    patch(...
        [Xmat(iii,2)+d/2,Xmat(iii,2)+d/2,Xmat(iii,2)-d/2,Xmat(iii,2)-d/2],...
        [Ymat(iii,2)-d/2,Ymat(iii,2)+d/2,Ymat(iii,2)+d/2,Ymat(iii,2)-d/2],...
        [Zmat(iii,2),Zmat(iii,2),Zmat(iii,2),Zmat(iii,2)],...
        colormat(colorind,:))
    colorind = colorind+1;
end

% Floor
for iii = 1:(n^2)
    patch(...
        [Xmat(iii,3)+d/2,Xmat(iii,3)+d/2,Xmat(iii,3)-d/2,Xmat(iii,3)-d/2],...
        [Ymat(iii,3)-d/2,Ymat(iii,3)+d/2,Ymat(iii,3)+d/2,Ymat(iii,3)-d/2],...
        [Zmat(iii,3),Zmat(iii,3),Zmat(iii,3),Zmat(iii,3)],...
        colormat(colorind,:))
    colorind = colorind+1;
end

% Right-hand-side wall
for iii = 1:(n^2)
    patch(...
        [Xmat(iii,4),Xmat(iii,4),Xmat(iii,4),Xmat(iii,4)],...
        [Ymat(iii,4)+d/2,Ymat(iii,4)+d/2,Ymat(iii,4)-d/2,Ymat(iii,4)-d/2],...
        [Zmat(iii,4)-d/2,Zmat(iii,4)+d/2,Zmat(iii,4)+d/2,Zmat(iii,4)-d/2],...
        colormat(colorind,:))
    colorind = colorind+1;
end

% Left-hand-side wall
for iii = 1:(n^2)
    patch(...
        [Xmat(iii,5),Xmat(iii,5),Xmat(iii,5),Xmat(iii,5)],...
        [Ymat(iii,5)+d/2,Ymat(iii,5)+d/2,Ymat(iii,5)-d/2,Ymat(iii,5)-d/2],...
        [Zmat(iii,5)-d/2,Zmat(iii,5)+d/2,Zmat(iii,5)+d/2,Zmat(iii,5)-d/2],...
        colormat(colorind,:))
    colorind = colorind+1;
end

% Camera settings
set(gca,'CameraPosition',[0 -5 0],'CameraTarget',[0 0 0],'CameraViewAngle',10)

% Axis settings
axis off

