% SH_LM_metrics_part2
%
% diffuseness metric by Xia et al., 2016
% brilliance metric by Fan et al., 2018
%
% The example1 of this code calculate the diffuseness metric and brilliance metric
% for high resolution light maps (Debevec, 1998) from USC database (http://gl.ict.usc.edu/data/highresprobes/)
% to select light maps that represent the ambient, focus, and brilliance lighting modes the best
%
% The example2 of this code tests the effects of contrast/gamma
% variation for "Eucalyptus" environment (Motoyoshi & Matoba 2012)
% 
% 2018(c) FZ
clc;
clear;
close all;


luminance_gamma = [1,0.5,0.06];
luminance_contrast = [0.95 1;0.75 1];

d = cell(5,1);
Img=hdrread('Eucalyptus_crop.hdr'); % to test effects of contrast/gamma variation in Motoyoshi 2012
luminance_map = 0.2126.*Img(:,:,1) + 0.7152.*Img(:,:,2) +0.0722.*Img(:,:,3);    % the conversion equations for linear signals
% figure,imshow(luminance_map),

for num=1:5 % number of maps being tested

    if num<=2
        I = imadjust(luminance_map,[],luminance_contrast(num,:),1);
    else
%       
        I = imadjust(luminance_map,[],[],1./luminance_gamma(num-2));
    end

%     figure(2),subplot(5,1,num),imshow(I), % uncomment this line to show light maps after gamma/contrast changes


    map_size=size(I);
    num_pixel = map_size(1)*map_size(2);
    

    %% initialize light source
    size1=map_size(1); %num of row
    size2=map_size(2); %num of coloum 
         
    luminance_matrix=zeros(size1,size2,3);
   
    for i=1:size1
        theta_s=i/size1*pi;
        
        for j=1:size2
            phi_s=(2*pi/size2)*(j-(size2+1)/2);
            
            luminance=I(i,j)*sin(theta_s);
            luminance_matrix(i,j,:)=[theta_s, phi_s, luminance];
        end
    end
        
    P=luminance_matrix(:,:,3);
    
    %% first 11 SH calculation (0th - 10th) 
    
    num_order = 11; % zero-th order included, thus "the highest order" = num-order -1
    
    % below FZ modified LX's code to calculate higher orders SH
    
    % Note: k.m code invloves command: factoral(N), which is only accurate under N<=21, and
    % here N <= l+abs(m) <= 2*l, so "order" should be no more than 10 in this
    % numerical method, i.e. use no more than num_order = 11 
    
     

      for order = 1:num_order
          l = order - 1;
          num_basis_func = 2 * l + 1;
          if l~=0  
          M{order} = [-l:1:l];
          else
          M{order} = 0;
          end
  
      
          for i=1:size(M{order},2)

              Ylm(:,:,i,order)=y_l_m(l,M{order}(i),luminance_matrix(:,:,1), luminance_matrix(:,:,2));
              f_value(:,:,i,order)=Ylm(:,:,i,order).*(P);
              Clm_square(i) = (sum(sum(f_value(:,:,i,order))))^2;
              d{num}(order) = sqrt(sum(Clm_square)); % SH_l, l = order - 1 the strength of each mode l
          end
          
      end
      w(num,1:3) = [sum(d{num}(1)) / sum(d{num}); sum(d{num}(2)) / sum(d{num}); sum(d{num}(4:num_order)) / sum(d{num})];
      w0 = w(num,1:3);
      f_maps(num) = sum(sum(sum(sum(f_value))));
end


  for i = 1:5
        D(i,:) = d{i};
  end
  
  
figure, 
hold on, Dx_Br=[];

for i = 1:5 % number of total maps
    Dx = 1-D(i,2)/D(i,1)/1.732; % Diffuseness metric (Xia et al., 2016)
    Br = sum(D(i,4:length(D)))./sum(D(i,:),2); % our brilliance metric: % 4:length(D) --> 4 here means starting with the third orders, skip clamp lighting (2nd order)
    Dx_Br = [Dx_Br;[Dx Br]]; 
end
    bar(Dx_Br');
legend('contrast 0.05','contrast 0.25','original','gamma 0.5','gamma 0.06'), grid on
set(gca, 'XTick',1:2,'XTickLabel',['Diffuseness metric (Xia et al. 2016)';'        our Brilliance metric       ' ]);grid off;
set(gca, 'FontSize',20);

