MATLAB Examples

function [lensV lensF]  = chordLengthDistribution(img,N,pxMin,calLen)

%Chord Length Distribution from Binary 2D Images

%Measure the chord length distribution of two phases in binary images

% A chord is a line segment with every point in one of the two phases,
% and with its end points on the interface between the two phases [1].
% Randomly oriented lines can be laid across each image, and the segments
% with uninterrupted intersections with either of two phases, p1 or p2,
% are tabulated. An implementation of this measurement based on
% the Bresenham's line algorithm was used for our paper [2].

%There are two considerations for this implementation:
% i.  Some directional anisotropy may be introduced by the method (randi)
%     used to generate random lines.
% ii. Technically, chords touching the edge of the images should be
%     removed; but I left them in to permit display of both phases from
%     the sample images.

% References:
% [1] S. Torquato, 1993, J. Chem. Phys. 98, 6472.
% [2] MacIver, M.R., Pawlik, M. 2017 Chem. Eng. Technol. 1–10.
%     DOI:10.1002/ceat.201600523

% INPUT
% img: a logical image (1's and 0's)
% N: number of line segments to sample (= < 10000)
% pxMin:minimum length of chords in pixels
% calLen: calibrated pixel centre spacing

% OUTPUT
% lensV: array of chords in phase 2 (p2, voids, 0's)
% lensF: array of chords in phase 1(p1, occupied, 1's)

% size of img
sz = size(img);

% random points in X and Y direction
X = randi(sz(1),[N,1,2]);
Y = randi(sz(2),[N,1,2]);
PT = [X Y];

% find LINES for BOX intersection
LINE = createLine(PT(:,:,1),PT(:,:,2));

% setup bounding box
XMIN = 1;
XMAX = sz(1);
YMIN = 1;
YMAX = sz(2);

BOX = [XMIN XMAX YMIN YMAX];

edge = clipLine(LINE, BOX);

% Sample image values along line segments
for a = 1:N
    [Xb{a} Yb{a}] = bresenham(edge(a,1),edge(a,2),edge(a,3),edge(a,4)); % get coordinates
    linearIndex{a} = sub2ind(sz,Xb{a},Yb{a}); % get index
    strp{a} = img(linearIndex{a}); % extract line
end

% analyze line segments and separate to agg and void chords
for aa = 1:N
    [lengths{aa}, values{aa}] = getChords(strp{aa}');
    lensF2{aa} = lengths{aa}(values{aa}==1);
    lensV2{aa} = lengths{aa}(values{aa}==0);
end

% collect into single array
if isempty(lensF2) == 0
    lensF = cat(2,lensF2{:});
elseif isempty(lensF2) == 1
    lensF = NaN;
end

if isempty(lensV2) == 0
    lensV = cat(2,lensV2{:});
elseif isempty(lensV2) == 1
    lensV = NaN;
end

% filter chords based on pxMin
lensV(lensV<=pxMin == 1) = [];
lensV(isnan(lensV) == 1)=[];
lensV = lensV*calLen;

lensF(isnan(lensF) == 1)=[];
lensF(lensF<=pxMin == 1) = [];
lensF = lensF*calLen;


Published with MATLAB® R2016a