Skip to content

Commit 0a7f466

Browse files
dleglanddlegland
committed
imGranulometry: add imOrientedGranulo.m for orientation granulometry
1 parent bc558fb commit 0a7f466

File tree

5 files changed

+217
-9
lines changed

5 files changed

+217
-9
lines changed

‎matImage/imGranulometry/Contents.m‎

Lines changed: 11 additions & 8 deletions
Original file line numberDiff line numberDiff line change
@@ -6,16 +6,19 @@
66
%
77
%
88
% Computation of granulometry curves
9-
% imGranulo - Compute granulometry curve of a given image.
10-
% imGranuloByRegion - Granulometry curve for each region of label image.
9+
% imGranulo - Compute granulometry curve of a given image.
10+
% imGranuloByRegion - Granulometry curve for each region of label image.
1111
%
1212
% Analysis of granulometry curves
13-
% granuloMeanSize - Compute geometric mean of granulometric curve.
14-
% granuloMean - Compute arithmetic mean of granulometric curve(s).
15-
% granuloStd - Compute standard deviation of granulometric curve(s).
16-
%
17-
% Variations of granulometry analysis
18-
% imOrientedGranulo - Gray level granulometry mean size for various orientations.
13+
% granuloMeanSize - Compute geometric mean of granulometric curve.
14+
% granuloMean - Compute arithmetic mean of granulometric curve(s).
15+
% granuloStd - Compute standard deviation of granulometric curve(s).
16+
%
17+
% Oriented granulometry analysis
18+
% imDirectionalGranulo - Directional granulometries for several orientations.
19+
% imGranuloOrientationMap - Orientation map of directional granulometry.
20+
% orientedLineStrel - Create an oriented line structuring element.
21+
% imOrientedGranulo - Gray level granulometry mean size for various orientations.
1922
%
2023
%
2124
% References

‎matImage/imGranulometry/imDirectionalGranulo.m‎

Lines changed: 95 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,95 @@
1+
function [res, orientList] = imDirectionalGranulo(img, nOrients, grType, LMax, varargin)
2+
% Directional granulometries for several orientations.
3+
%
4+
% Usage:
5+
% RES = imDirectionalGranulo(IMG, NORIENT, GRTYPE, LMAX)
6+
%
7+
% Input parameters:
8+
% IMG: input image, 2D gray-level
9+
% NORIENT: the number of orientations to consider
10+
% GRTYPE: the type of granulomtry (only 'opening' is supported for now)
11+
% LMAX: the maximum size of line
12+
%
13+
%
14+
% Example
15+
% imDirectionalGranulo
16+
%
17+
% See also
18+
% orientedLineStrel
19+
%
20+
% Reference
21+
% The methodology is described in the following article:
22+
% "Oriented granulometry to quantify fibre orientation distributions in
23+
% synthetic and plant fibre composite preforms", by V. Gager, D. Legland,
24+
% A. Bourmaud, A. Le Duigou, F. Pierre, K. Behlouli and C. Baley. (2020),
25+
% Industrial Crops and Products 152, p. 112548.
26+
% doi: https://doi.org/10.1016/j.indcrop.2020.112548
27+
%
28+
29+
30+
%
31+
32+
% ------
33+
% Author: David Legland
34+
% e-mail: [email protected]
35+
% Created: 2018-12-18, using Matlab 9.5.0.944444 (R2018b)
36+
% Copyright 2018 INRA - Cepia Software Platform.
37+
38+
39+
%% Input arguments
40+
41+
if ~strcmp(grType, 'opening')
42+
error('Only ''opening'' granulometry type is currently supported');
43+
end
44+
45+
46+
%% Initialization
47+
48+
dim = size(img);
49+
50+
% create the list of angles
51+
orientList = linspace(0, 180, nOrients+1);
52+
orientList(end) = [];
53+
54+
% create the list of line diameters
55+
% (consider only odd values to ensure symetry of the strel)
56+
diamList = 1:2:LMax;
57+
nSteps = length(diamList);
58+
59+
% allocate memory for global result
60+
res = zeros([dim nOrients], 'double');
61+
62+
% allocate memory for intermediate results
63+
resOp = zeros([dim nSteps+1], 'double');
64+
65+
% image for normalizing granulometry curves
66+
refImage = double(img);
67+
refImage(img <= 0) = 1;
68+
69+
70+
%% Main processing
71+
72+
% iterate over orientations
73+
for iOrient = 1:nOrients
74+
disp(sprintf('Orient: %d/%d', iOrient, nOrients)); %#ok<DSPS>
75+
76+
% angles from horizontal, in degrees, in CW order
77+
% (correspond to CCW when visualizing image results)
78+
theta = -orientList(iOrient);
79+
80+
% iterate over granulometry steps to create a stack of results
81+
for i = 1:nSteps
82+
se = orientedLineStrel(diamList(i), theta);
83+
resOp(:,:,i) = imopen(img, se);
84+
end
85+
86+
% compute granulometry curve for each pixel
87+
gr = bsxfun(@rdivide, cat(3, zeros(dim), diff(-resOp, 1, 3)), refImage) * 100;
88+
89+
% compute mean size for each position
90+
meanSizes = granuloMeanSize(reshape(gr, [numel(img) nSteps+1]), [diamList LMax]);
91+
92+
% stores the mean size for the current orientation
93+
res(:,:,iOrient) = reshape(meanSizes, size(img));
94+
end
95+

‎matImage/imGranulometry/imGranuloOrientationMap.m‎

Lines changed: 44 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,44 @@
1+
function [orientMap, resMax, rgb] = imGranuloOrientationMap(img, nOrients, grType, LMax)
2+
% Orientation map of directional granulometry.
3+
%
4+
% ORIMAP = imGranuloOrientationMap(IMG, NORIENT, GRTYPE, LMAX)
5+
%
6+
% [ORIMAP, RGB] = imGranuloOrientationMap(IMG, NORIENT, GRTYPE, LMAX)
7+
% Also returns an RGB version for display
8+
%
9+
% Example
10+
% imGranuloOrientationMap
11+
%
12+
% See also
13+
% imGranulometry, imDirectionalGranulo, granuloMeanSize
14+
%
15+
16+
% ------
17+
% Author: David Legland
18+
% e-mail: [email protected]
19+
% Created: 2018-12-19, using Matlab 9.5.0.944444 (R2018b)
20+
% Copyright 2018 INRA - Cepia Software Platform.
21+
22+
% compute the 3D array of mean size for each position and each orientation
23+
[res, orientList] = imDirectionalGranulo(img, nOrients, grType, LMax);
24+
resMax = max(res, [], 3);
25+
26+
% convert to radians, over all directions
27+
orientRads = deg2rad(2 * orientList);
28+
29+
% compute average direction for each pixel, weighted by granulometric size
30+
dxList = reshape(cos(orientRads), [1 1 nOrients]);
31+
dyList = reshape(sin(orientRads), [1 1 nOrients]);
32+
dxMoy = sum(bsxfun(@times, res, dxList), 3) ./ sum(res, 3);
33+
dyMoy = sum(bsxfun(@times, res, dyList), 3) ./ sum(res, 3);
34+
35+
% create orientation map, in degrees
36+
orientMap = mod(rad2deg(atan2(dyMoy, dxMoy) / 2) + 180, 180);
37+
38+
% optionnaly create an rgb version
39+
if nargout > 1
40+
hue = orientMap / 180;
41+
sat = double(img) / double(max(img(:)));
42+
val = ones(size(img));
43+
rgb = hsv2rgb(cat(3, hue, sat, val));
44+
end

‎matImage/imGranulometry/imOrientedGranulo.m‎

Lines changed: 2 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -19,7 +19,8 @@
1919
% imOrientedGranulo
2020
%
2121
% See also
22-
% imGranulometry, imGranulo, imGranuloByRegion, granuloMeanSize
22+
% imGranulometry, imDirectionalGranulo, imGranulo, imGranuloByRegion,
23+
% granuloMeanSize
2324

2425
% ------
2526
% Author: David Legland

‎matImage/imGranulometry/orientedLineStrel.m‎

Lines changed: 65 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,65 @@
1+
function se = orientedLineStrel(L, theta)
2+
% Create an oriented line structuring element.
3+
%
4+
% SE = orientedLineStrel(L, THETA)
5+
% Generates a binary images corresponding the linear structuring element
6+
% with length L and orientation THETA (in degrees).
7+
% The length corresponds to the approwimated Euclidean length of the
8+
% final structuring element.
9+
%
10+
% Example
11+
% % Creates a structuring element with length 10 pixels and 30 degrees
12+
% % orientation
13+
% se = orientedLineStrel(10, 30)
14+
% se =
15+
% 5x9 logical array
16+
%
17+
% 1 1 0 0 0 0 0 0 0
18+
% 0 0 1 1 0 0 0 0 0
19+
% 0 0 0 0 1 0 0 0 0
20+
% 0 0 0 0 0 1 1 0 0
21+
% 0 0 0 0 0 0 0 1 1
22+
%
23+
% See also
24+
% imGranulometry, imDirectionalGranulo, imGranulo
25+
%
26+
27+
% ------
28+
% Author: David Legland
29+
% e-mail: [email protected]
30+
% Created: 2018-12-18, using Matlab 9.5.0.944444 (R2018b)
31+
% Copyright 2018 INRA - Cepia Software Platform.
32+
33+
% pre-compute trigonometric quantities
34+
cost = cos(deg2rad(theta));
35+
sint = sin(deg2rad(theta));
36+
37+
% compute strel size
38+
if abs(cost) >= abs(sint)
39+
% horizontal strel directions
40+
xRadius = round((abs(L * cost) - 1) / 2);
41+
yRadius = round(xRadius * abs(sint / cost));
42+
else
43+
% vertical strel directions
44+
yRadius = round((abs(L * sint) - 1) / 2);
45+
xRadius = round(yRadius * abs(cost / sint));
46+
end
47+
48+
% allocate memory
49+
dim = [2*yRadius+1 2*xRadius+1];
50+
se = false(dim);
51+
52+
if abs(cost) >= abs(sint)
53+
% Process horizontal strel directions
54+
lx = -xRadius:xRadius;
55+
ly = round(lx * sint / cost);
56+
inds = sub2ind(dim, ly + yRadius + 1, lx + xRadius + 1);
57+
se(inds) = true;
58+
59+
else
60+
% Process vertical strel directions
61+
ly = -yRadius:yRadius;
62+
lx = round(ly * cost / sint);
63+
inds = sub2ind(dim, ly + yRadius + 1, lx + xRadius + 1);
64+
se(inds) = true;
65+
end

0 commit comments

Comments
 (0)