diff --git a/matImage/imGranulometry/Contents.m b/matImage/imGranulometry/Contents.m index 2e80154..3b05fe8 100644 --- a/matImage/imGranulometry/Contents.m +++ b/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 diff --git a/matImage/imGranulometry/imDirectionalGranulo.m b/matImage/imGranulometry/imDirectionalGranulo.m new file mode 100644 index 0000000..65ca680 --- /dev/null +++ b/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: david.legland@inra.fr +% 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 + + % 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 + diff --git a/matImage/imGranulometry/imGranuloOrientationMap.m b/matImage/imGranulometry/imGranuloOrientationMap.m new file mode 100644 index 0000000..7ad84b4 --- /dev/null +++ b/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: david.legland@inrae.fr +% 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 diff --git a/matImage/imGranulometry/imOrientedGranulo.m b/matImage/imGranulometry/imOrientedGranulo.m index d622cd4..a74b5bd 100644 --- a/matImage/imGranulometry/imOrientedGranulo.m +++ b/matImage/imGranulometry/imOrientedGranulo.m @@ -19,7 +19,8 @@ % imOrientedGranulo % % See also -% imGranulometry, imGranulo, imGranuloByRegion, granuloMeanSize +% imGranulometry, imDirectionalGranulo, imGranulo, imGranuloByRegion, +% granuloMeanSize % ------ % Author: David Legland diff --git a/matImage/imGranulometry/orientedLineStrel.m b/matImage/imGranulometry/orientedLineStrel.m new file mode 100644 index 0000000..1ab9b03 --- /dev/null +++ b/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: david.legland@inrae.fr +% 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