Skip to content

Instantly share code, notes, and snippets.

@zabela

zabela/getGlobalContrastFactor.m

Created January 21, 2014 12:29
Show Gist options
  • Star 1 You must be signed in to star a gist
  • Fork 2 You must be signed in to fork a gist
  • Save zabela/8539136 to your computer and use it in GitHub Desktop.
Save zabela/8539136 to your computer and use it in GitHub Desktop.
Download ZIP
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
Raw
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
@Arun-Niranjan
Copy link

At line 34, is it an issue that the resolutions variable isn't used? So in actual fact the superpixel resolutions used by this function are [1 2 4 8 16 32 64 128 256]?

Perhaps this should be changed to

rIm = imresize(im, 1/(resolutions(i)-1), 'bilinear');

But I may have mis-understood?

@rayryeng
Copy link

rayryeng commented Jun 6, 2017

@Arun-Niranjan I agree. It should be changed to how you specified it. I just read the original paper.

@rayryeng
Copy link

rayryeng commented Jun 7, 2017 •
edited

This code is unfortunately fraught with bugs. There are if/else conditions that don't make any sense. For example, if (x == 1) && (x == H) is never satisfied because x can never be both the first and last row at the same time. There is also the case where the resolution array is not used. The code below fixes this as well as removes the unnecessary if statements. I've also added a comment on what LC is doing for the output in the comments block:

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 
%       LC - local contrast at each scale
%

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

@rayryeng
Copy link

rayryeng commented Jun 7, 2017 •
edited

Another implementation that is more vectorized. This is better for older versions of MATLAB, but after R2016b, the for loop approach and the vectorized implementations are nearly equal in performance. Take note that there will be some slight computational differences. The old code resized the image originally in unsigned integer prior to normalizing and finding the perceptual luminance. This method converts to double first for precision reasons, then does the resizing after.

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');
    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

@dvolgyes
Copy link

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 ?

@vik748
Copy link

vik748 commented May 7, 2020

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

@Jdissopa
Copy link

thank you for this code

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment