GLCM for texture analysis • GLCM (gray level co-occurrence matrix) is mainly useful to perform the texture analysis and find the features from image.
• As name suggested its work on gray image and try to create sort of 2 d histogram from image.
• Main application of GLCM are texture analysis, feature extraction and segmentation.

Steps to Calculate GLCM matrix :

• Lets assume image I which is gray image.
• Initialize GLCM matrix size 256 x 256  (256 is level of GLCM).
• Suppose we use zero angle of GLCM means direction of GLCM is horizontal.
• Suppose distance of GLCM is 1, means we just look horizontally next pixel to current pixel.
• In image  I(i,j) get the gray value (suppose value of pixel is a =  127 at I(i,j)), and get gray  value I(i,j+1) in case of distance 1 and zero degree of  GLCM (suppose value of pixel is b = 58 at I(i,j+1)).
•  Go to GLCM matrix co-ordinate (a = 127, b=58) and increment the value by 1.
• Iterate the full image that will give us GLCM matrix of zero degree for distance 1.
• According to texture type, GLCM distance can be decided.
• For road  because texture changes so rapidly we consider the small distance to calculate GLCM but for bricks, distance of the GLCM matrix should be large.
• Angle of GLCM should be selected according to  direction of image texture changes, so in brick image we want to consider 2 direction 0 degree and 90 degree.
• Calculate features from GLCM matrix, there are many features are available to perform the texture analysis like contrast, correlation, energy , and homogeneity etc.
• We will get the contrast feature from GLCM matrix which are sufficient to say weather texture is rough or smooth.

Lets write basic GLCM code which calculate the zero degree GLCM for 256 level and get contrast feature.

%%%%%%%%%%%%%%%%%%%%%%%%%Code Start here %%%%%%%%%%%%%%%

%%%%%  GLCM function start here %%%%%%%%%%%%%

function GLCM_0 = getGLCM0(img_gray, distance)
%% function calculate the GLCM matrix at zero degree angle and given distance
%% GLCM_0 is GLCM matrix for 0 degree angle.
%% img_gray is input gray image
%% distance is distance of GLCM calculated

%% initialize GLCM matrix

GLCM_0 = zeros(256,256);

[row,col] = size(img_gray);

for i=1:row
for j=1:col-distance

pixel1 = img_gray(i,j)+1;
pixel2 = img_gray(i,j+distance)+1;
GLCM_0(pixel1,pixel2) = GLCM_0(pixel1,pixel2)+1;

end
end

end

%%% GLCM function end here %%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%% Contrast function start here %%%%%%%%%%%%%%%%%

function Contrast = GLCM_contrast(GLCM)

%% function calculate the GLCM contrast for GLCM matrix
%% GLCM is Gray level co-occurance matrix as input.
%% Contasrt is feature of GLCM matrix

%% normalize the GLCM

[row,col] = size(GLCM);

N_factor = sum(sum(GLCM));

GLCM = GLCM./N_factor;

Contrast = 0;

for i = 1:row
for j = 1:col

temp = ((i-j)^2)* GLCM(i,j);
Contrast = Contrast + temp;
end
end

end

%%%% Contrast function end here  %%%%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%main GLCM code %%%%%%%%%

%% this main code calls  GLCM and contrast function for given image

clc % clear the screen

clear all % clear all variables

close all % close all MATLAB windows

% set path to read image from system

ImagePath = ‘D:\DSUsers\uidp6927\image_processingCode\frontend-large Brick.jpg’;

img_gray = rgb2gray(img_RGB); % Convert image to gray

img_gray = imresize(img_gray,[900,1200]);   % re sizing (normalizing the size)

GLCM_0 = getGLCM0(img_gray, 5);     %% get the GLCM matrix for distance 5

Contrast = GLCM_contrast(GLCM_0);   %% calculate the contrast feature

disp(‘Contrast of function is =’)
Contrast
%% use matlab inbuilt function to verify result

glcm = graycomatrix(img_gray,’NumLevels’,256,’offset’,[0 5]);

stats = graycoprops(GLCM_0,’Contrast’);

disp(‘Contrast of MATLAB is = ‘)
stats.Contrast

%%%%%%%%%%%%%%%%%%%%%%%%Code End here   %%%%%%%%%%%%%%%%%

Matlab provide the inbuilt function from which we can calculate the GLCM and Contrast and Verify our results.

glcm = graycomatrix(img_gray,’NumLevels’,256,’offset’,[0 5]);

stats = graycoprops(GLCM_0,’Contrast’);

Lest take two input image first is road texture and second is brick image apply the GLCM and get the contrast feature.

Result:

Contrast feature of figure. 2 road texture is  3.9489e+03

Contrast feature of figure.3 Brick texture is 91.0206

As from result we can say that road texture have so much texture because of that its contrast value is so much high if less texture will be present in image contrast value will be low.

So for smooth area in image contrast will be low,  for rough area contrast will be high.

* Note: GLCM value depend on the size of image we are performing the GLCM so size of image should be fix or normalize.

For more  information about GLCM  feature analysis please go through the paper harlick texture features .  Happy Learning

Cheers

2 D Projective Transform Lets take the Lena image to perform the transformation

• In this post I am going to discuss image projective transforms.
• These transform are easy to understand and useful in many applications like object detection,  3 D view geometry and ADAS.
• Suppose image I(x,y) as input image and our output (transform) image is O(u,v)
•  So we can write a Transform matrix T such that :

O(u,v) = T [I(x,y)]

• Suppose image is in homogeneous co-ordinate system we can add 1

O(u,v,1) = T [I(x,y,1)]

• Where T is 3 x 3 matrix, according to T matrix we can able to transform  the image many ways.

Lets write a code which can take input image multiply it with transform matrix and produce transform image

%%%%%%%%%%%%%%%%%%%%%%%% Code Start%%%%%%%%%%%%%%%%

function transform_img = trasform_2d(T, input_img)

%function perform the 2 D projective transform

% input_image : gray input image

% T is 3 x 3 trasnform matrix

%transform_image : output gray transform image

[row,col] = size(input_img);

transform_img = zeros(2*row,2*col);

for i=-row:2*row
for j=-col:2*col

U = inv(T)*[i j 1]’;      % where ‘ is transpose

if(U(3)>1)                       % in case of 8 dof projective transform

U(1) = U(1)/U(3); %% taking care homogeneous co ordinate
U(2) = U(2)/U(3);
end

if((U(1)>0)&&(U(1)<row-1)&&(U(2)>0)&&(U(2)<col-1)) %% consider only points                                                                                                                 %which lie inside of input image

transform_img(round(i+row+1),round(j+col+1)) =                                         input_img(round(U(1)+1),round(U(2)+1));

end

end
end

end

%%%%%%%%%%%% Code End %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

Following are the main 2 D project transform :

1.  Euclidean Transform
2. Similarity Transform
3. Affine Transform
4. Projective Transform

Euclidean Transform :

Let define transformation matrix for Euclidean transform

T = [ R11  R12  Tx  ;  R21  R22  Ty ; 0 0 1]

Where R11 , R12, R21 and  R22 are rotation matrix and Tx and Ty are translation. So it have have 3 degree of freedom 1 in rotation and 2 in translation.

lets take 30 degree rotation angle and -100 pixel translation in both direction Result of Euclidean Transform

Euclidean Transform dont have property to scale the images.

Similarity Transform:

Let define transformation matrix for Similarity transform

T = [ s*R11  s*R12  Tx  ; s* R21  s*R22  Ty ; 0 0 1]

It have 4 degree of freedom.

if in Euclidean transform we added the feature of scale (s) its become similarity tranform

lets take 40 degree rotation angle and -100 pixel translation in both direction and 1.5 scaling factor

Affine Transform :

Let define transformation matrix for Affine transform

T = [ a11  a12  Tx  ; a21  a22  Ty ; 0 0 1]

In Similarity transform if we can able to take different angle as well as different scale it will become affine transform.

lets take 40 degree rotation angle in x direction and 30 degree rotation angle in y direction , -100 pixel translation in both direction, 0.7 scaling factor in x direction and 1.5  scaling factor in y direction.

Projective Transform :

Define the transformation full matrix which include the homogeneous  factor as well.

T =  [h11  h12  h13 ; h21  h22  h23 ; h31  h32  h33]

Projective transform have 8 degree of freedom because we added 2  coefficient h31  h32 in affine transform matrix which take care of image homogeneous co-ordinate geometry. (might be in furure calibration post i will discuss it about in detail)

%%%%%%%%%%%%%%%%%%  Code Start  %%%%%%%%%%%%%%%%%%%%%%%

clc % clear the screen

clear all % clear all variables

close all % close all MATLAB windows

% set path to read image from system

ImagePath = ‘D:\DSUsers\uidp6927\image_processingCode\lena.jpg’;

img_RGB = im2double(img_RGB); % Convert image to double precision

img_gray = rgb2gray(img_RGB); % convert image to gray

[row, col] = size(img_gray);

figure,imshow(img_gray)
%For a nonreflective similarity,
%[u,v] = [x,y,1]T

%% euclidean transformation.

scale = 1; % scale factor
angle = 30*pi/180; % rotation angle
tx = -100; % x translation
ty = -100; % y translation

sc = scale*cos(angle);
ss = scale*sin(angle);

T = [ sc ss tx;
-ss sc ty;
0 0 1];
% call the function transform_2d
transform_img = trasform_2d(T, img_gray);
figure, imshow(transform_img)
title(‘Euclidean transform’)

%% Similarity transform
scale = 1.5; % scale factor
angle = 40*pi/180; % rotation angle
tx = -100; % x translation
ty = -50; % y translation

sc = scale*cos(angle);
ss = scale*sin(angle);

T = [ sc ss tx;
-ss sc ty;
0 0 1];

transform_img = trasform_2d(T, img_gray);
figure, imshow(transform_img)
title(‘Similarity transform’)
%% Affine transform
scale1 = .7; % scale factor
scale2 = 1.7;
angle1 = 40*pi/180; % rotation angle
angle2 = 30*pi/180; % rotation angle
tx = -100; % x translation
ty = -50; % y translation
sc1 = scale1*cos(angle1);
ss1 = scale1*sin(angle2);
sc2 = scale2*cos(angle1);
ss2 = scale2*sin(angle2);

T = [ sc1 ss1 tx;
-ss2 sc2 ty;
0 0 1];
transform_img = trasform_2d(T, img_gray);

figure, imshow(transform_img)
title(‘Affine transform’)

%% projective transform
scale1 = .7; % scale factor
scale2 = 1.7;
angle1 = 40*pi/180; % rotation angle
angle2 = 30*pi/180; % rotation angle
tx = -100; % x translation
ty = -50; % y translation

sc1 = scale1*cos(angle1);
ss1 = scale1*sin(angle2);
sc2 = scale2*cos(angle1);
ss2 = scale2*sin(angle2);

T = [ sc1 ss1 tx;
-ss2 sc2 ty;
.01 .02 1];

transform_img = trasform_2d(T, img_gray);
figure, imshow(transform_img)
title(‘Projective transform’)

%%%%%%%%%%%%%%%%%%%%%  Code End %%%%%%%%%%%%%%%%%%%%%%%  Happy Learning

Cheers

Wavelet Transform In this post, I will discuss the Haar wavelet analysis (wavelet transform) in image processing.

Applications of Haar wavelet transform in image processing are Feature extraction, Texture analysis, Image compression etc.

• Wavelet is time frequency resolutions analysis.
• As shown in above figure.1  its have the information of different time frequencies of signal.
• Its like applying a comb of filters bands in signal they are called levels. First we apply full frequency band and have no knowledge of time.
• Apply filter which have half bandwidth of signal lose some frequency information but  have some knowledge of time resolution called first level decomposition.
• Further reduce the filter bandwidth will be second level decomposition.

Haar wavelet is made of two filters one is low pass filter and second high pass filter.

Coefficients of Low pass filter  = [0.707, 0.707]

Coefficients of High pass filter  = [-0.707, 0.707]

• Figure .3 shows the block diagram of first level decomposition.
• I is image of size m , n. We take the image column vectors (1d) and apply the low pass filter and high pass filter.
•  Down sample the output and reshape the image I1 of size m, n/2 (size is n/2 because we down sampled the columns by 2).
• Apply low and high pass filter in row wise in image I1 again down sample the output and reshape the image (m/2,n/2).

• Figure. 5 shows the result of haar wavelet decomposition results at first level.
• LL filter output will proceeds for next level decomposition. As its contains most information of input image and looks almost same at down sampled resolution.
• HH, LH  and HL will have texture information in vertical, horizontal and ramp directions. HH, LH and HL are called haar features.
• We can further  decompose the LL image into LL1, LH1, HL1 and HH1 images.

MATLAB provide Wavelet Toolbox GUI to perform the wavelet analysis.

Type wevemenu in command window of  MATLAB.

%%%%%%Code Start %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function [LL, LH, HL, HH] = haar_wave(image_in)

%% function haar_wave calculate the haar wave decomposition of given imput image
%% Image in :- input gray level image
%% LL :- Low Low band output image
%% LH :- Low High band output image
%% HL :- High Low band output image
%% HH :- High High band output image

image_in = im2double(image_in);
filter_L = [.707, .707];
filter_H = [-.707, .707];
[row, col] = size(image_in);
image_vector = image_in(:);

%% apply low pass filter in one d image column wise image
image_low_pass_1D = filter(filter_L,1,image_vector);
image_high_pass_1D = filter(filter_H,1,image_vector);

% Down sample the image low and high pass image.
temp = 1;
Dimage_low_pass_1D = 0;
Dimage_high_pass_1D = 0;

for i = 1:2:length(image_low_pass_1D)
Dimage_low_pass_1D(temp) = image_low_pass_1D(i);
Dimage_high_pass_1D(temp) = image_high_pass_1D(i);
temp = temp+1;
end

%% reshape the image in 2 D, size will be row/2 , col
Dimage_low_pass = reshape(Dimage_low_pass_1D,round(row/2), col)’;
Dimage_high_pass = reshape(Dimage_high_pass_1D, round(row/2),col)’;

%% conver the image in 1 D but in column vector, I transpose Dimage_low_paas image above.

Dimage_low_pass_1D = Dimage_low_pass(:);
Dimage_high_pass_1D = Dimage_high_pass(:);

LL_1D = filter(filter_L,1,Dimage_low_pass_1D);
LH_1D = filter(filter_H,1,Dimage_low_pass_1D);

HL_1D = filter(filter_L,1,Dimage_high_pass_1D);
HH_1D = filter(filter_H,1,Dimage_high_pass_1D);

%Down sample the all 1 _d vector images
temp = 1;
DLL_1D = 0;
DLH_1D = 0;
DHL_1D = 0;
DHH_1D = 0;

for i = 1:2:length(LL_1D)
DLL_1D(temp) = LL_1D(i);
DLH_1D(temp) = LH_1D(i);
DHL_1D(temp) = HL_1D(i);
DHH_1D(temp) = HH_1D(i);
temp = temp+1;
end

%% reshape the images in size row/2 and clo/2

LL = reshape(DLL_1D, round(col/2), round(row/2))’;
LH = reshape(DLH_1D, round(col/2), round(row/2))’;
HL = reshape(DHL_1D, round(col/2), round(row/2))’;
HH = reshape(DHH_1D, round(col/2), round(row/2))’;

end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%% call the haar wavelet function and check the result

clc % clear the screen

clear all % clear all variables

close all % close all MATLAB windows

% set path to read image from system

ImagePath = ‘D:\DSUsers\uidp6927\image_processingCode\lena.jpg’;

img_RGB = im2double(img_RGB); % Convert image to double precision

img_gray = rgb2gray(img_RGB); % convert image to gray

% first level decomposition
[LL, LH, HL, HH] = haar_wave(img_gray);

% second level decomposition
[LL1, LH1, HL1, HH1] = haar_wave(LL);

%% Draw all image

figure, imshow(img_gray);

figure,
subplot(2,2,1), imshow(LL,[])
title(‘LL haar wave image’)
subplot(2,2,2), imshow(LH,[])
title(‘LH haar wave image’)
subplot(2,2,3), imshow(HL,[])
title(‘HL haar wave image’)
subplot(2,2,4), imshow(HH,[])
title(‘HH haar wave image’)

figure,
subplot(2,2,1), imshow(LL1,[])
title(‘LL1 haar wave image’)
subplot(2,2,2), imshow(LH1,[])
title(‘LH1 haar wave image’)
subplot(2,2,3), imshow(HL1,[])
title(‘HL1 haar wave image’)
subplot(2,2,4), imshow(HH1,[])
title(‘HH1 haar wave image’)

%%%%%%%%%%%%%%%%%Code End %%%%%%%%%%%%%%%%%%%%%%%%%%%%  Happy Learning

Cheers