imGranulometry: add imOrientedGranulo.m for orientation granulometry

matImage/imGranulometry/Contents.m

@@ -6,16 +6,19 @@
 %
 %
 % Computation of granulometry curves
-%   imGranulo         - Compute granulometry curve of a given image.
-%   imGranuloByRegion - Granulometry curve for each region of label image.
+%   imGranulo               - Compute granulometry curve of a given image.
+%   imGranuloByRegion       - Granulometry curve for each region of label image.
 %
 % Analysis of granulometry curves
-%   granuloMeanSize   - Compute geometric mean of granulometric curve.
-%   granuloMean       - Compute arithmetic mean of granulometric curve(s).
-%   granuloStd        - Compute standard deviation of granulometric curve(s).
-%
-% Variations of granulometry analysis
-%   imOrientedGranulo - Gray level granulometry mean size for various orientations.
+%   granuloMeanSize         - Compute geometric mean of granulometric curve.
+%   granuloMean             - Compute arithmetic mean of granulometric curve(s).
+%   granuloStd              - Compute standard deviation of granulometric curve(s).
+%
+% Oriented granulometry analysis
+%   imDirectionalGranulo    - Directional granulometries for several orientations.
+%   imGranuloOrientationMap - Orientation map of directional granulometry.
+%   orientedLineStrel       - Create an oriented line structuring element.
+%   imOrientedGranulo       - Gray level granulometry mean size for various orientations.
 %
 %
 % References

matImage/imGranulometry/imDirectionalGranulo.m

@@ -0,0 +1,95 @@
+function [res, orientList] = imDirectionalGranulo(img, nOrients, grType, LMax, varargin)
+% Directional granulometries for several orientations.
+%
+%   Usage:
+%   RES = imDirectionalGranulo(IMG, NORIENT, GRTYPE, LMAX)
+%
+%   Input parameters:
+%   IMG:         input image, 2D gray-level
+%   NORIENT:     the number of orientations to consider
+%   GRTYPE:      the type of granulomtry (only 'opening' is supported for now)
+%   LMAX:        the maximum size of line
+%
+%
+%   Example
+%     imDirectionalGranulo
+%
+%   See also
+%     orientedLineStrel
+%
+%   Reference
+%   The methodology is described in the following article:
+%   "Oriented granulometry to quantify fibre orientation distributions in
+%   synthetic and plant fibre composite preforms", by V. Gager, D. Legland,
+%   A. Bourmaud, A. Le Duigou, F. Pierre, K. Behlouli and C. Baley. (2020),
+%   Industrial Crops and Products 152, p. 112548. 
+%   doi: https://doi.org/10.1016/j.indcrop.2020.112548
+%
+
+
+%
+
+% ------
+% Author: David Legland
+% e-mail: [email protected]
+% Created: 2018-12-18,    using Matlab 9.5.0.944444 (R2018b)
+% Copyright 2018 INRA - Cepia Software Platform.
+
+
+%% Input arguments
+
+if ~strcmp(grType, 'opening')
+    error('Only ''opening'' granulometry type is currently supported');
+end
+
+
+%% Initialization
+
+dim = size(img);
+
+% create the list of angles
+orientList = linspace(0, 180, nOrients+1);
+orientList(end) = [];
+
+% create the list of line diameters 
+% (consider only odd values to ensure symetry of the strel)
+diamList = 1:2:LMax;
+nSteps = length(diamList);
+
+% allocate memory for global result
+res = zeros([dim nOrients], 'double');
+
+% allocate memory for intermediate results
+resOp = zeros([dim nSteps+1], 'double');
+
+% image for normalizing granulometry curves
+refImage = double(img);
+refImage(img <= 0) = 1;
+
+
+%% Main processing 
+
+% iterate over orientations
+for iOrient = 1:nOrients
+    disp(sprintf('Orient: %d/%d', iOrient, nOrients)); %#ok<DSPS>
+
+    % angles from horizontal, in degrees, in CW order 
+    % (correspond to CCW when visualizing image results)
+    theta = -orientList(iOrient);
+
+    % iterate over granulometry steps to create a stack of results
+    for i = 1:nSteps
+        se = orientedLineStrel(diamList(i), theta);
+        resOp(:,:,i) = imopen(img, se);
+    end
+
+    % compute granulometry curve for each pixel
+    gr = bsxfun(@rdivide, cat(3, zeros(dim), diff(-resOp, 1, 3)), refImage) * 100;
+
+    % compute mean size for each position
+    meanSizes = granuloMeanSize(reshape(gr, [numel(img) nSteps+1]), [diamList LMax]);
+
+    % stores the mean size for the current orientation
+    res(:,:,iOrient) = reshape(meanSizes, size(img));
+end
+

matImage/imGranulometry/imGranuloOrientationMap.m

@@ -0,0 +1,44 @@
+function [orientMap, resMax, rgb] = imGranuloOrientationMap(img, nOrients, grType, LMax)
+% Orientation map of directional granulometry.
+%
+%   ORIMAP = imGranuloOrientationMap(IMG, NORIENT, GRTYPE, LMAX)
+%
+%   [ORIMAP, RGB] = imGranuloOrientationMap(IMG, NORIENT, GRTYPE, LMAX)
+%   Also returns an RGB version for display
+%
+%   Example
+%   imGranuloOrientationMap
+%
+%   See also
+%     imGranulometry, imDirectionalGranulo, granuloMeanSize
+%
+
+% ------
+% Author: David Legland
+% e-mail: [email protected]
+% Created: 2018-12-19,    using Matlab 9.5.0.944444 (R2018b)
+% Copyright 2018 INRA - Cepia Software Platform.
+
+% compute the 3D array of mean size for each position and each orientation
+[res, orientList] = imDirectionalGranulo(img, nOrients, grType, LMax);
+resMax = max(res, [], 3);
+
+% convert to radians, over all directions
+orientRads = deg2rad(2 * orientList);
+
+% compute average direction for each pixel, weighted by granulometric size
+dxList = reshape(cos(orientRads), [1 1 nOrients]);
+dyList = reshape(sin(orientRads), [1 1 nOrients]);
+dxMoy = sum(bsxfun(@times, res, dxList), 3) ./ sum(res, 3);
+dyMoy = sum(bsxfun(@times, res, dyList), 3) ./ sum(res, 3);
+
+% create orientation map, in degrees
+orientMap = mod(rad2deg(atan2(dyMoy, dxMoy) / 2) + 180, 180);
+
+% optionnaly create an rgb version
+if nargout > 1
+    hue = orientMap / 180;
+    sat = double(img) / double(max(img(:)));
+    val = ones(size(img));
+    rgb = hsv2rgb(cat(3, hue, sat, val));
+end

matImage/imGranulometry/imOrientedGranulo.m

@@ -19,7 +19,8 @@
 %   imOrientedGranulo
 %
 %   See also
-%     imGranulometry, imGranulo, imGranuloByRegion, granuloMeanSize
+%     imGranulometry, imDirectionalGranulo, imGranulo, imGranuloByRegion,
+%     granuloMeanSize 
 
 % ------
 % Author: David Legland

matImage/imGranulometry/orientedLineStrel.m

@@ -0,0 +1,65 @@
+function se = orientedLineStrel(L, theta)
+% Create an oriented line structuring element.
+%
+%   SE = orientedLineStrel(L, THETA)
+%   Generates a binary images corresponding the linear structuring element
+%   with length L and orientation THETA (in degrees).
+%   The length corresponds to the approwimated Euclidean length of the
+%   final structuring element.
+%
+%   Example
+%     % Creates a structuring element with length 10 pixels and 30 degrees
+%     % orientation
+%     se = orientedLineStrel(10, 30)
+%     se =
+%       5x9 logical array
+% 
+%        1   1   0   0   0   0   0   0   0
+%        0   0   1   1   0   0   0   0   0
+%        0   0   0   0   1   0   0   0   0
+%        0   0   0   0   0   1   1   0   0
+%        0   0   0   0   0   0   0   1   1
+%
+%   See also
+%     imGranulometry, imDirectionalGranulo, imGranulo
+%
+
+% ------
+% Author: David Legland
+% e-mail: [email protected]
+% Created: 2018-12-18,    using Matlab 9.5.0.944444 (R2018b)
+% Copyright 2018 INRA - Cepia Software Platform.
+
+% pre-compute trigonometric quantities
+cost = cos(deg2rad(theta));
+sint = sin(deg2rad(theta));
+
+% compute strel size
+if abs(cost) >= abs(sint)
+    % horizontal strel directions
+    xRadius = round((abs(L * cost) - 1) / 2);
+    yRadius = round(xRadius * abs(sint / cost));
+else
+    % vertical strel directions
+    yRadius = round((abs(L * sint) - 1) / 2);
+    xRadius = round(yRadius * abs(cost / sint));
+end
+
+% allocate memory
+dim = [2*yRadius+1 2*xRadius+1];
+se = false(dim);
+
+if abs(cost) >= abs(sint)
+    % Process horizontal strel directions
+    lx = -xRadius:xRadius;
+    ly = round(lx * sint / cost);
+    inds = sub2ind(dim, ly + yRadius + 1, lx + xRadius + 1);
+    se(inds) = true;
+
+else
+    % Process vertical strel directions
+    ly = -yRadius:yRadius;
+    lx = round(ly * cost / sint);
+    inds = sub2ind(dim, ly + yRadius + 1, lx + xRadius + 1);
+    se(inds) = true;
+end