% SH_LM_metrics_example1
%
% 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
%
% Also see SH_LM_metrics_example2. it tests the effects of contrast/gamma variation for "Eucalyptus" environment (Motoyoshi & Matoba 2012)
% 
% 2018(c) FZ
clc;
clear;
close all;

% Imagename=cell(1,13);
% Imagename={'pisa.hdr','glacier.hdr','doge2.hdr','uffizi-large.hdr','grace-new.hdr','ennis.hdr','Park.hdr','Parking lot.hdr','River_Road.hdr','Woods.hdr','City_Night_Lights.hdr','Stonewall.hdr','Tunnel.hdr','night_free.hdr','Container.hdr','Harbor.hdr'}
num_maps=6;
Imagename=cell(1,num_maps);
Imagename={'glacier.hdr','pisa.hdr','doge.hdr','grace-new.hdr','uffizi-large.hdr','ennis.hdr'};
d = cell(num_maps,1);
luminance_map_all = zeros(1,num_maps);

% figure,imshow(luminance_map),

for map=1:num_maps % number of maps
    
    Imageused=Imagename{map};
    Img=hdrread(Imageused);

    luminance_map = 0.2126.*Img(:,:,1) + 0.7152.*Img(:,:,2) +0.0722.*Img(:,:,3);    % the conversion equations for linear signals (Stokes et al., 2012)

%     I = luminance_map;
%     figure(2),subplot(num_maps,1,map),imshow(luminance_map), % show all light maps

    map_size=size(luminance_map);
    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=luminance_map(i,j)*sin(theta_s);
            luminance_matrix(i,j,:)=[theta_s, phi_s, luminance];
        end
    end
        
    P=luminance_matrix(:,:,3);
    
    %% first x orders SH calculation (0th - (x-1)th) 
    
    num_order = 11; % zero-th order included, thus "the highest order" = num-order -1
    
    % num_order =  2 is sufficient for diffuseness metric D_xia
    
    % 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))).*(1/2500))^2; 
              % here i removed .*(1/2500) in the original code from Xia, it probably has something to do with the size of the map she tested using the original code
              % removing this won't change the metrics score
              Clm_square(i) = (sum(sum(f_value(:,:,i,order))))^2; %
              d{map}(order) = sqrt(sum(Clm_square)); % SH_l, l = order - 1 the strength of each mode l
          end
          
      end
      w(map,1:3) = [sum(d{map}(1)) / sum(d{map}); sum(d{map}(2)) / sum(d{map}); sum(d{map}(4:num_order)) / sum(d{map})];
      w0 = w(map,1:3);
      f_maps(map) = sum(sum(sum(sum(f_value))));
end


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

for i = 1:num_maps % number of total maps
    Dx = 1-D(i,2)/D(i,1)/sqrt(3); % 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('pisa.hdr','glacier.hdr','doge.hdr','uffizi-large.hdr','grace-new.hdr','ennis.hdr'), grid on
set(gca, 'XTick',1:2,'XTickLabel',['Diffuseness metric (Xia et al. 2016)  ';'Brilliance metric (Zhang et al., 2018)' ]);grid off;
set(gca, 'FontSize',20);

