-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathglrlm.m
144 lines (117 loc) · 4.27 KB
/
glrlm.m
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
function Final_Out = glrlm(img)
% [SRE,LRE,GLN,RP,RLN,LGRE,HGRE]= glrlm(img,16,mask)
% gray level run length matrix computation
%
% input:
% - img, an input grayscale image (RGB images are converted to grayscale)
% - quantize, quantization levels. Normally set to 16. Should be larger than 1.
% - mask, a binary mask to use with values of 1 at the ROI's. If whole image
% is needed then mask = ones(size(img(:,:,1)));
%
% output: texture features
%
% 1. SHORT RUN EMPHASIS (SRE)
% 2. LONG RUN EMPHASIS(LRE)
% 3. GRAY LEVEL NON-UNIFORMITY (GLN)
% 4. RUN PERCENTAGE (RP)
% 5. RUN LENGTH NON-UNIFORMITY (RLN)
% 6. LOW GRAY LEVEL RUN EMPHASIS (LGRE)
% 7. HIGH GRAY LEVEL RUN EMPHASIS (HGRE)
%
% example:
% I = imread('cameraman.tif');
% imshow(I)
% mask = ones(size(I(:,:,1)));
% quantize = 16;
% [SRE,LRE,GLN,RP,RLN,LGRE,HGRE] = glrlm(I,quantize,mask)
mask = ones(size(img(:,:,1)));
quantize = 16;
% if color => make gray scale
if size(img,3)>1
img = im2gray(img);
end
img = im2double(img); % to double
% crop the image to the mask bounds for faster processing
stats = regionprops(mask,'BoundingBox');
bx = int16(floor(stats.BoundingBox)) + int16(floor(stats.BoundingBox)==0);
img = img(bx(2):bx(2)+bx(4)-1,bx(1):bx(1)+bx(3)-1);
mask = mask(bx(2):bx(2)+bx(4)-1,bx(1):bx(1)+bx(3)-1);
% adjust range
mini = min(img(:)); % find minimum
img = img-mini; % let the range start at 0
maxi = max(img(:)); % find maximum
% quantize the image to discrete integer values in the range 1:quantize
levels = maxi/quantize:maxi/quantize:maxi-maxi/quantize;
img = imquantize(img,levels);
% apply the mask
img(~mask) = 0;
% initialize glrlm: p(i,j)
% - with i the amount of bin values (quantization levels)
% - with j the maximum run length (because yet unknown, assume maximum length
% of image)
% - four different orientations are used (0, 45, 90 and 135 degrees)
p0 = zeros(quantize,max(size(img)));
p45 = zeros(quantize,max(size(img)));
p90 = zeros(quantize,max(size(img)));
p135 = zeros(quantize,max(size(img)));
% initialize maximum value for j
maximgS = max(size(img));
% add zeros to the borders
img = padarray(img,[1 1]);
% initialize rotation
img45 = imrotate(img,45);
% find the run length for each quantization level
for i = 1:quantize
% find the pixels corresponding to the quantization level
BW = int8(img == i);
BWr = int8(img45 == i);
% find the start and end points of the run length
G0e = (BW(2:end-1,2:end-1) - BW(2:end-1,3:end)) == 1;
G0s = (BW(2:end-1,2:end-1) - BW(2:end-1,1:end-2)) == 1;
G45e = (BWr(2:end-1,2:end-1) - BWr(2:end-1,3:end)) == 1;
G45s = (BWr(2:end-1,2:end-1) - BWr(2:end-1,1:end-2)) == 1;
G90e = (BW(2:end-1,2:end-1) - BW(3:end,2:end-1)) == 1;
G90s = (BW(2:end-1,2:end-1) - BW(1:end-2,2:end-1)) == 1;
G135e = (BWr(2:end-1,2:end-1) - BWr(3:end,2:end-1)) == 1;
G135s = (BWr(2:end-1,2:end-1) - BWr(1:end-2,2:end-1)) == 1;
% find the indexes
G0s = G0s'; G0s = find(G0s(:));
G0e = G0e'; G0e = find(G0e(:));
G45s = G45s'; G45s = find(G45s(:));
G45e = G45e'; G45e = find(G45e(:));
G90s = find(G90s(:));
G90e = find(G90e(:));
G135s = find(G135s(:));
G135e = find(G135e(:));
% find the lengths
lengths0 = G0e - G0s + 1;
lengths45 = G45e - G45s + 1;
lengths90 = G90e - G90s + 1;
lengths135 = G135e - G135s + 1;
% fill the matrix
p0(i,:) = hist(lengths0,1:maximgS);
p45(i,:) = hist(lengths45,1:maximgS);
p90(i,:) = hist(lengths90,1:maximgS);
p135(i,:) = hist(lengths135,1:maximgS);
end
% add all orientations
p = p0 + p45 + p90 + p135;
% calculate the features
Out0 = ComputePar(p0,maximgS,quantize,mask);
Out45 = ComputePar(p45,maximgS,quantize,mask);
Out90 = ComputePar(p90,maximgS,quantize,mask);
Out135 = ComputePar(p135,maximgS,quantize,mask);
Out = ComputePar(p,maximgS,quantize,mask);
Final_Out = [Out0; Out45; Out90; Out135; Out];
end
function out = ComputePar(p,maximgS,quantize,mask)
totSum = sum(p(:));
SRE = sum(sum(p,1) ./ ((1:maximgS).^2)) / totSum;
LRE = sum(sum(p,1) .* ((1:maximgS).^2)) / totSum;
RLN = sum(sum(p,1) .^2) / totSum;
RP = totSum / sum(mask(:));
GLN = sum(sum(p,2) .^2) / totSum;
LGRE = sum(sum(p,2) .* ((1:quantize)'.^2)) / totSum;
HGRE = sum(sum(p,2) .^2) / totSum;
out = [SRE,LRE,GLN,RP,RLN,LGRE,HGRE]';
end