% radialcenter.m
%
% Copyright 2011-2012, Raghuveer Parthasarathy, The University of Oregon
%
%%
% Disclaimer / License
% This program is free software: you can redistribute it and/or
% modify it under the terms of the GNU General Public License as
% published by the Free Software Foundation, either version 3 of the
% License, or (at your option) any later version.
% This set of programs is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
% General Public License for more details.
% You should have received a copy of the GNU General Public License
% (gpl.txt) along with this program. If not, see .
%%
%
% Calculates the center of a 2D intensity distribution.
% Method: Considers lines passing through each half-pixel point with slope
% parallel to the gradient of the intensity at that point. Considers the
% distance of closest approach between these lines and the coordinate
% origin, and determines (analytically) the origin that minimizes the
% weighted sum of these distances-squared.
%
% Inputs
% I : 2D intensity distribution (i.e. a grayscale image)
% Size need not be an odd number of pixels along each dimension
%
% Outputs
% xc, yc : the center of radial symmetry,
% px, from px #1 = left/topmost pixel
% So a shape centered in the middle of a 2*N+1 x 2*N+1
% square (e.g. from make2Dgaussian.m with x0=y0=0) will return
% a center value at x0=y0=N+1.
% Note that y increases with increasing row number (i.e. "downward")
% sigma : Rough measure of the width of the distribution (sqrt. of the
% second moment of I - min(I));
% Not determined by the fit -- output mainly for consistency of
% formatting compared to my other fitting functions, and to get
% an estimate of the particle "width"
%
% Raghuveer Parthasarathy
% The University of Oregon
% August 21, 2011 (begun)
% last modified Apr. 6, 2012 (minor change)
% Copyright 2011-2012, Raghuveer Parthasarathy
function [xc yc sigma] = radialcenter(I)
% Number of grid points
[Ny Nx] = size(I);
% grid coordinates are -n:n, where Nx (or Ny) = 2*n+1
% grid midpoint coordinates are -n+0.5:n-0.5;
% The two lines below replace
% xm = repmat(-(Nx-1)/2.0+0.5:(Nx-1)/2.0-0.5,Ny-1,1);
% and are faster (by a factor of >15 !)
% -- the idea is taken from the repmat source code
xm_onerow = -(Nx-1)/2.0+0.5:(Nx-1)/2.0-0.5;
xm = xm_onerow(ones(Ny-1, 1), :);
% similarly replacing
% ym = repmat((-(Ny-1)/2.0+0.5:(Ny-1)/2.0-0.5)', 1, Nx-1);
ym_onecol = (-(Ny-1)/2.0+0.5:(Ny-1)/2.0-0.5)'; % Note that y increases "downward"
ym = ym_onecol(:,ones(Nx-1,1));
% Calculate derivatives along 45-degree shifted coordinates (u and v)
% Note that y increases "downward" (increasing row number) -- we'll deal
% with this when calculating "m" below.
dIdu = I(1:Ny-1,2:Nx)-I(2:Ny,1:Nx-1);
dIdv = I(1:Ny-1,1:Nx-1)-I(2:Ny,2:Nx);
% Smoothing --
h = ones(3)/9; % simple 3x3 averaging filter
fdu = conv2(dIdu, h, 'same');
fdv = conv2(dIdv, h, 'same');
dImag2 = fdu.*fdu + fdv.*fdv; % gradient magnitude, squared
% Slope of the gradient . Note that we need a 45 degree rotation of
% the u,v components to express the slope in the x-y coordinate system.
% The negative sign "flips" the array to account for y increasing
% "downward"
m = -(fdv + fdu) ./ (fdu-fdv);
% *Very* rarely, m might be NaN if (fdv + fdu) and (fdv - fdu) are both
% zero. In this case, replace with the un-smoothed gradient.
NNanm = sum(isnan(m(:)));
if NNanm > 0
unsmoothm = (dIdv + dIdu) ./ (dIdu-dIdv);
m(isnan(m))=unsmoothm(isnan(m));
end
% If it's still NaN, replace with zero. (Very unlikely.)
NNanm = sum(isnan(m(:)));
if NNanm > 0
m(isnan(m))=0;
end
%
% Almost as rarely, an element of m can be infinite if the smoothed u and v
% derivatives are identical. To avoid NaNs later, replace these with some
% large number -- 10x the largest non-infinite slope. The sign of the
% infinity doesn't matter
try
m(isinf(m))=10*max(m(~isinf(m)));
catch
% if this fails, it's because all the elements are infinite. Replace
% with the unsmoothed derivative. There's probably a more elegant way
% to do this.
m = (dIdv + dIdu) ./ (dIdu-dIdv);
end
% Shorthand "b", which also happens to be the
% y intercept of the line of slope m that goes through each grid midpoint
b = ym - m.*xm;
% Weighting: weight by square of gradient magnitude and inverse
% distance to gradient intensity centroid.
sdI2 = sum(dImag2(:));
xcentroid = sum(sum(dImag2.*xm))/sdI2;
ycentroid = sum(sum(dImag2.*ym))/sdI2;
w = dImag2./sqrt((xm-xcentroid).*(xm-xcentroid)+(ym-ycentroid).*(ym-ycentroid));
% least-squares minimization to determine the translated coordinate
% system origin (xc, yc) such that lines y = mx+b have
% the minimal total distance^2 to the origin:
% See function lsradialcenterfit (below)
[xc yc] = lsradialcenterfit(m, b, w);
%%
% Return output relative to upper left coordinate
xc = xc + (Nx+1)/2.0;
yc = yc + (Ny+1)/2.0;
% A rough measure of the particle width.
% Not at all connected to center determination, but may be useful for tracking applications;
% could eliminate for (very slightly) greater speed
Isub = I - min(I(:));
[px,py] = meshgrid(1:Nx,1:Ny);
xoffset = px - xc;
yoffset = py - yc;
r2 = xoffset.*xoffset + yoffset.*yoffset;
sigma = sqrt(sum(sum(Isub.*r2))/sum(Isub(:)))/2; % second moment is 2*Gaussian width
%%
function [xc yc] = lsradialcenterfit(m, b, w)
% least squares solution to determine the radial symmetry center
% inputs m, b, w are defined on a grid
% w are the weights for each point
wm2p1 = w./(m.*m+1);
sw = sum(sum(wm2p1));
smmw = sum(sum(m.*m.*wm2p1));
smw = sum(sum(m.*wm2p1));
smbw = sum(sum(m.*b.*wm2p1));
sbw = sum(sum(b.*wm2p1));
det = smw*smw - smmw*sw;
xc = (smbw*sw - smw*sbw)/det; % relative to image center
yc = (smbw*smw - smmw*sbw)/det; % relative to image center
end
end