Algorithm to measure image contrast, adapted from "Global contrast factor-a new approach to image contrast" (Matkovic, Kresimir et al., 2005) http://citeseerx.ist.psu.edu/viewdoc/download?doi=10.1.1.84.2683&rep=rep1&type=pdf

getGlobalContrastFactor.m

function [GCF, LC] = getGlobalContrastFactor( im )
%
% GCF = getGlobalContrastFactor( im )
%   
% MATLAB algorithm implementation of the 
% "Global contrast factor-a new approach to image contrast" 
% (Matkovic, Kresimir et al., 2005) 
%
% http://www.cg.tuwien.ac.at/research/publications/2005/matkovic-2005-glo/
%
%   Input:
%       im  - image in grayscale
%
%   Output:
%       GCF - global contrast factor 
%

    % 9 different resolution levels
    GCF = 0.0;
    
    resolutions = [1 2 4 8 16 25 50 100 200];
    
    LC = zeros(size(resolutions));
    
    W = size(im,2);
    H = size(im,1);
    
    rIm = im;
    
    for i=1:length(resolutions)
        
        %attempt at resizing as in the paper
        if i>1
          rIm = imresize(im, 1/(2^(i-1)), 'bilinear');      
        end

        W = size(rIm,2);
        H = size(rIm,1);
        
        rL = zeros(size(rIm));
        % compute linear luminance l
        l = (double(rIm(:,:))/255) * 2.2;

        % compute perceptual luminance L
        rL(:,:) = 100 * sqrt(l);
        
        % compute local contrast for each pixel
        lc = 0.0;
        for x=1:H
            for y=1:W
                
                if (x == 1) && (x == H)
                    if (y == 1) && (y == W)
                       lc = lc + 0;
                    elseif (y == 1)
                       lc = lc + abs(rL(x, y) - rL(x,y+1));
                    elseif (y == W)
                        lc = lc + abs(rL(x, y) - rL(x,y-1));
                    else
                        lc = lc + ( abs(rL(x, y) - rL(x,y-1)) + ...
                                abs(rL(x, y) - rL(x,y+1)) )/2;
                    end
                                
                elseif (x == 1)
                    if (y == 1) && (y == W) 
                        lc = lc + abs(rL(x, y) - rL(x+1,y));
                    elseif (y == 1) 
                        lc = lc + ( abs(rL(x, y) - rL(x,y+1)) + ...
                                abs(rL(x, y) - rL(x+1,y)) )/2;
                    elseif (y == W) 
                        lc = lc + ( abs(rL(x, y) - rL(x,y-1)) + ...
                                abs(rL(x, y) - rL(x+1,y)) )/2;
                    else 
                        lc = lc + ( abs(rL(x, y) - rL(x,y-1)) + ...
                                abs(rL(x, y) - rL(x,y+1)) + ...
                                abs(rL(x, y) - rL(x+1,y)) )/3;
                    end
                    
                elseif (x == H)
                    if (y == 1) && (y == W)
                        lc = lc + abs(rL(x, y) - rL(x-1,y));
                    elseif (y == 1) 
                        lc = lc + ( abs(rL(x, y) - rL(x,y+1)) + ...
                                abs(rL(x, y) - rL(x-1,y)) )/2;
                    elseif (y == W) 
                        lc = lc + ( abs(rL(x, y) - rL(x,y-1)) + ...
                                abs(rL(x, y) - rL(x-1,y)) )/2;
                    else 
                        lc = lc + ( abs(rL(x, y) - rL(x,y-1)) + ...
                                    abs(rL(x, y) - rL(x,y+1)) + ...
                                    abs(rL(x, y) - rL(x-1,y)) )/3;
                    end
                else % x > 1 && x < H
                    if (y == 1) && (y == W)
                        lc = lc + ( abs(rL(x, y) - rL(x+1,y)) + ...
                                abs(rL(x, y) - rL(x-1,y)) )/2;
                    elseif (y == 1) 
                        lc = lc + ( abs(rL(x, y) - rL(x,y+1)) + ...
                                    abs(rL(x, y) - rL(x+1,y)) + ...
                                    abs(rL(x, y) - rL(x-1,y)) )/3;
                    elseif (y == W) 
                        lc = lc + ( abs(rL(x, y) - rL(x,y-1)) + ...
                                    abs(rL(x, y) - rL(x+1,y)) + ...
                                    abs(rL(x, y) - rL(x-1,y)) )/3;
                    else 
                        lc = lc + ( abs(rL(x, y) - rL(x,y-1)) + ...
                                    abs(rL(x, y) - rL(x,y+1)) + ...
                                    abs(rL(x, y) - rL(x-1,y)) + ...
                                    abs(rL(x, y) - rL(x+1,y)) )/4;
                    end

                end
            end
        end
        
        % compute average local contrast c
        c(i) = lc/(W*H);
        w(i) = (-0.406385*(i/9)+0.334573)*(i/9)+ 0.0877526;
        
        % compute global contrast factor
        LC(i) = c(i)*w(i);
        GCF = GCF + LC(i);
        
    end  
end

Maybe I miss something, but following the original article, some parts seem to be strange, e.g.

   % Compute perceptual luminance
    l_resize = 100 * ((im_resize * 2.2).^(0.5));

How is the gamma correction performed? Instead of *2.2, shouldn't it be ^2.2 ?

Hi, thanks for putting this out there. Matlab's imresize does antialiasing by default leading to slightly different results from what is described in the paper. Also taking into account what is pointed out by @dvolgyes here is another version of the code:

function out = globalContrastFactor(im)

superpixels = [1 2 4 8 16 25 50 100 200];

% Convert colour to grayscale
if size(im, 3) ~= 1
    im = rgb2gray(im);
end

% Convert to double precision and normalize
im = im2double(im);

% Stores resized image
im_resize = im;

out = 0; % Output GCF measure

% For each superpixel scale...
for ii = 1 : 9
    
    % Resize the image as per the paper
    if ii ~= 1
        im_resize = imresize(im, 1.0 / superpixels(ii), 'bilinear','Antialiasing',false);
    end
    
    % Compute perceptual luminance
    l_resize = 100 * ((im_resize .^ 2.2).^(0.5));

    % Determine an array of scales where when we sum up the local
    % differences, we divide by a certain value.  4 is when we have
    % all valid neighbours, 3 is when we're at a border and 2 is when
    % we are at any of the four corners.
    scale = 4 * ones(size(l_resize));
    scale(:, [1 end]) = 3;
    scale([1 end], :) = 3;    
    scale([1 end], [1 end]) = 2;
    
    % Pad the image with a one pixel border
    l_pad = zeros(size(l_resize) + 2);
    l_pad(2:end-1,2:end-1) = l_resize;
    
    % Calculate the local differences
    left = abs(l_resize - l_pad(2:end-1,1:end-2));
    up = abs(l_resize - l_pad(1:end-2,2:end-1));
    right = abs(l_resize - l_pad(2:end-1, 3:end));
    down = abs(l_resize - l_pad(3:end, 2:end-1));
    
    % Zero out those locations that do not have valid neighbours
    left(:, 1) = 0;
    up(1, :) = 0;
    right(:, end) = 0;
    down(end, :) = 0;
    
    % Calculate the local contrast for this superpixel scale    
    lc_scale = (left + right + up + down) ./ scale; 
    
    % Calculate the weight at this superpixel scale
    wi = (-0.406385 * (ii / 9) + 0.334573) * (ii / 9) + 0.0877526;
    
    % Now calculate the GCF at this     
    out = out + wi * mean(lc_scale(:));
end
end

Also here is an equivalent version in python:

import cv2
import numpy as np

def compute_image_average_contrast(k, gamma=2.2):
    L = 100 * np.sqrt((k / 255) ** gamma )
    # pad image with border replicating edge values
    L_pad = np.pad(L,1,mode='edge')

    # compute differences in all directions
    left_diff = L - L_pad[1:-1,:-2]
    right_diff = L - L_pad[1:-1,2:]
    up_diff = L - L_pad[:-2,1:-1]
    down_diff = L - L_pad[2:,1:-1]

    # create matrix with number of valid values 2 in corners, 3 along edges and 4 in the center
    num_valid_vals = 3 * np.ones_like(L)
    num_valid_vals[[0,0,-1,-1],[0,-1,0,-1]] = 2
    num_valid_vals[1:-1,1:-1] = 4

    pixel_avgs = (np.abs(left_diff) + np.abs(right_diff) + np.abs(up_diff) + np.abs(down_diff)) / num_valid_vals

    return np.mean(pixel_avgs)

def compute_global_contrast_factor(img):
    if img.ndim != 2:
        gr = cv2.cvtColor(img, cv2.COLOR_BGR2GRAY)
    else:
        gr = img

    superpixel_sizes = [1, 2, 4, 8, 16, 25, 50, 100, 200]

    gcf = 0

    for i,size in enumerate(superpixel_sizes,1):
        wi =(-0.406385 * i / 9 + 0.334573) * i/9 + 0.0877526
        im_scale = cv2.resize(gr, (0,0), fx=1/size, fy=1/size,
                              interpolation=cv2.INTER_LINEAR)
        avg_contrast_scale = compute_image_average_contrast(im_scale)
        gcf += wi * avg_contrast_scale

    return gcf

thank you for this code