-
-
Save zabela/8539136 to your computer and use it in GitHub Desktop.
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 I agree. It should be changed to how you specified it. I just read the original paper.
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
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
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
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?