Code covered by the BSD License  

Highlights from
Muscle fascicle tracking - Ultrasound

image thumbnail

Muscle fascicle tracking - Ultrasound



02 Sep 2011 (Updated )

Implementation of an optical flow algorithm to track muscle length changes imaged with ultrasound.

gsmooth2(im, sigmas, varargin)
function [im, reg] = gsmooth2(im, sigmas, varargin)
%GSMOOTH2 Gaussian image smoothing
%   [SMOOTH, REGION] = GSMOOTH2(IMAGE, SIGMAS) carries out Gaussian
%   smoothing on an image.
%   SIGMAS specifies the smoothing constants. This may be a matrix of the
%   form [SIGMAX SIGMAY] or a scalar specifying the values of both SIGMAX
%   and SIGMAY. SIGMAX and SIGMAY are the "sigma" parameters of the 1D
%   Gaussian masks used for smoothing along the rows and columns
%   respectively. (Smoothing along a row means smoothing across columns -
%   i.e. the mask is a row vector.)
%       The output arrays will be smaller than the input arrays in order to
%       avoid having to extrapolate beyond the boundaries of the array. The
%       result REGION reports the region of the input arrays for which the
%       output values have been estimated, in the form [MINROW, MAXROW,
%       MINCOL, MAXCOL]. The size of the output array is [MAXROW-MINROW+1,
%       MAXCOL-MINCOL+1]. The reduction in size depends on the smoothing
%       parameters, and is chosen so that the smoothing masks are a good
%       approximation to the untruncated Gaussian mask.
%   [SMOOTH, REGION] = GSMOOTH2(IMAGE, SIGMAS, WRAP) can be used to specify
%   that the image wraps round on one or more axes. WRAP may be a logical
%   1x2 matrix of the form [WRAPX WRAPY] or a logical scalar which
%   specifies both WRAPX and WRAPY. If WRAPX is true the rows of the images
%   wrap round; if WRAPY is true the columns wrap round. If the
%   rows/columns wrap, the number of columns/rows will not be reduced in
%   the output array. An empty WRAP is the same as omitting it, which is
%   the same as setting it to false.
%   region of interest to be specified. REGION may be a 4-element row
%   vector with elements [MINROW, MAXROW, MINCOL, MAXCOL] describing a
%   rectangular region. The results arrays will have size [MAXROW-MINROW+1,
%   MAXCOL-MINCOL+1] and will contain the estimated gradients for the
%   specified region of the input images. The REGION argument is returned
%   unchanged. An empty REGION is the same as omitting it.
%       It is possible to specify a region which goes right up to the
%       boundaries of the image, or even goes outside it. In these cases
%       reflection at the boundaries will be used to extrapolate the image
%       in any directions in which it does not wrap round.
%   [SMOOTH, REGION] = GSMOOTH2(IMAGE, SIGMAS, WRAP, 'same') specifies that
%   the output array should be the same size as the input array. Boundary
%   reflection is used to extend the input array prior to smoothing.
%   REGION = GSMOOTH2(IMAGE, SIGMAS, WRAP, 'region') returns only the
%   default region, without computing the gradients.

% Copyright David Young 2010

[sigmas, bcons, reg, convreg, ronly] = checkinputs(im, sigmas, varargin{:});

if ronly
    im = reg;
    if ~isempty(convreg)
        im = exindex(im, convreg(1):convreg(2), bcons{1}, ...
            convreg(3):convreg(4), bcons{2});
    im = gsmooth1(im, 1, sigmas(1));
    im = gsmooth1(im, 2, sigmas(2));


% -------------------------------------------------------------------------

function [sigmas, bcons, reg, convreg, regonly] = ...
    checkinputs(im, sigmas, wraps, reg)
% Check arguments and get defaults, plus input/output convolution regions
% and boundary conditions

error(nargchk(2, 4, nargin, 'struct'));

% image argument
validateattributes(im, {'double'}, {'2d'});

% sigmas argument
validateattributes(sigmas, {'double'}, {'nonnegative'});
if isscalar(sigmas)
    sigmas = [sigmas sigmas];
elseif isequal(size(sigmas), [1 2])
    sigmas = sigmas([2 1]);   % xy -> row col
    error('gsmooth:badsigmas', 'SIGMAS wrong size');

% wraps argument
if nargin < 3 || isempty(wraps)
    wraps = [false false];
    validateattributes(wraps, {'logical'}, {});
    if isscalar(wraps)
        wraps = [wraps wraps];
    elseif isequal(size(wraps), [1 2])
        wraps = wraps([2 1]); % xy -> row col for exindex
        error('gsmooth:badwraps', 'WRAPS wrong size');
boundopts = {'symmetric' 'circular'};
bcons = boundopts(wraps+1);

% region argument
if nargin < 4
    reg = [];
regonly = strcmp(reg, 'region');

imreg = [1 size(im,1) 1 size(im,2)];   % whole image region
mrg = gausshsize(sigmas);  % convolution margins
mrg = [1 -1 1 -1] .* mrg([1 1 2 2]);
if isempty(reg) || regonly
    % default region - small enough not to need extrapolation - shrink on
    % non-wrapped dimensions
    reg = imreg + ~wraps([1 1 2 2]) .* mrg;
elseif strcmp(reg, 'same')
    reg = imreg;
    validateattributes(reg, {'double'}, {'real', 'integer', 'size', [1 4]});
if any(reg([2 4]) < reg([1 3]))
    error('gsmooth:badreg', 'REGION or array size too small');
% compute input region for convolution - expand on all dimensions
convreg = reg - mrg;    % expand
if isequal(convreg, imreg)
    convreg = [];   % signal no trimming or padding


% -------------------------------------------------------------------------

function im = gsmooth1(im, dim, sigma)
% Smooth an image IM along dimension DIM with a 1D Gaussian mask of
% parameter SIGMA

hsize = gausshsize(sigma);  % reasonable truncation

msize = [1 1];
msize(dim) = 2*hsize+1;

if sigma > 0
    mask = fspecial('gauss', msize, sigma);
    im = conv2(im, mask, 'valid');


% -------------------------------------------------------------------------

function hsize = gausshsize(sigma)
% Default for the limit on a Gaussian mask of parameter sigma.
% Produces a reasonable degree of truncation without too much error.
hsize = ceil(2.6*sigma);

Contact us