% borrowed the idea from Steve's blog
% and wrote the function I needed, here is a link to Steve's blog
url = 'http://blogs.mathworks.com/images/steve/2010/rice_binary.png';
bw = imread(url);
dataArray = zeros(size(bw));
dataArray = rand(size(bw));
s = regionprops(bw, 'Orientation', 'MajorAxisLength', ...
'MinorAxisLength', 'Eccentricity', 'Centroid');
imagesc(bw)
hold on
phi = [0,pi/2, pi, 2*pi*3/4.0];
cosphi = cos(phi);
sinphi = sin(phi);
for k = 1:length(s)
xbar = s(k).Centroid(1);
ybar = s(k).Centroid(2);
a = s(k).MajorAxisLength/2;
b = s(k).MinorAxisLength/2;
theta = pi*s(k).Orientation/180;
R = [ cos(theta) sin(theta)
-sin(theta) cos(theta)];
xy = [a*cosphi; b*sinphi];
xy = R*xy;
x = xy(1,:) + xbar;
y = xy(2,:) + ybar;
plot(x,y,'*g');
[avg, totalElem, jx, jy] = avgInEllipse(xbar, ybar, a, b,...
s(k).Orientation, bw, dataArray);
xvec = zeros(jx.size(),1);
yvec = zeros(jy.size(),1);
for tt = 0:jx.size()-1
xvec(tt+1) = jx.elementAt(tt);
yvec(tt+1) = jy.elementAt(tt);
end
plot(xvec, yvec, '.m')
disp(['total no of interior points :' num2str(totalElem)]);
disp(['Average :' num2str(avg)]);
end
% this part is directly copied from Steve's blog, I just inserted comments
% create 50 points around the circumference of the unit circle
phi = linspace(0,2*pi,50);
cosphi = cos(phi);
sinphi = sin(phi);
for k = 1:length(s)
xbar = s(k).Centroid(1);
ybar = s(k).Centroid(2);
% get radii of the two circles centered at (0,0)
a = s(k).MajorAxisLength/2;
b = s(k).MinorAxisLength/2;
% convert Orientation angle to radian value
theta = pi*s(k).Orientation/180;
% rotation matrix
R = [ cos(theta) sin(theta)
-sin(theta) cos(theta)];
% a*cosphi gives points on the x-axis of a circle
% drawn around (0,0) point with a radius of 'a'
% b*sinphi gives points on the y-axis of a circle
% drawn around (0,0) point with a radius of 'b'
xy = [a*cosphi; b*sinphi];
% rotate the points to match orientation
xy = R*xy;
% translate the points that were generated around (0,0) point
% to the new center (xbar, ybar) aka centroid
x = xy(1,:) + xbar;
y = xy(2,:) + ybar;
plot(x,y,'r','LineWidth',2);
end
hold off
 |
| Input Image |
 |
| output image, enclosed points drawn in purple |
 |
output image, two elements are zoomed in
|
function [avg, totalElem, jvector_x, jvector_y] = avgInEllipse(cx, cy, a, b,...
orientation, bw, dataArray)
% Author : Mohammad A. Khan
% ======================= inputs ==========================
% this function is to be used in conjunction with regionprops
% cx, cy are the x and y coordinated of the centroid
% a, b are the 1/2 lengths of the major and minor axes of the
% ellipse returned by regionprops, orientation is the orientation
% returned by regionprops, which is the degree of angle between
% ellipse's major axis and x-axis, bw is the binary array where
% elements were detected, it can also be though of the mask array
% against the dataArray, this dataArray is the actual scalar
% array on which the averaging over ellipsoidal area is expected
% ==================== outputs =============================
% avg is the average over ellipse area calculated from dataArray
% totalElem is the total number of pixels in ellipse
% jvector_x, jvector_y are java Vectors of x,y coordinates
% of all the interior points of the ellipse
jvector_x = java.util.Vector;
jvector_y = java.util.Vector;
phi = [0,pi/2, pi, 2*pi*3/4.0];
cosphi = cos(phi);
sinphi = sin(phi);
theta = pi* orientation/180;
R = [ cos(theta) sin(theta)
-sin(theta) cos(theta)];
% calculate the four end points of major and minor axes
xy = [a*cosphi; b*sinphi];
% generate points on major axis
xx = ceil(min(xy(1,:))):floor(max(xy(1,:)));
totalElem = 0;
sumTotal = 0;
% for each point on major axis calculate all possible
% y locations that fall within ellipse
for i = xx
proot = [-sqrt((1-i*i/a/a)*b*b), sqrt((1-i*i/a/a)*b*b)];
yy = real((min(proot))):real((max(proot)));
xcopy = ones(1,length(yy)).*i;
% rotate the ellipse that was generated around (0,0) point
% to match the oreintation given and translate the ellipse to the
% centroid location
rtxy = R*[xcopy;yy];
xcor = (rtxy(1,:)+cx);
ycor = (rtxy(2,:)+cy);
% establish bound check
xout1 = find(xcor < 1);
xout2 = find(xcor > size(bw,2));
xcor([xout1, xout2]) = [];
ycor([xout1, xout2]) = [];
yout1 = find(ycor < 1);
yout2 = find(ycor > size(bw,1));
xcor([yout1, yout2]) = [];
ycor([yout1, yout2]) = [];
xcor = round(xcor);
ycor = round(ycor);
% delete all the points that fall on the background
% of the binary mask image
if ~isempty(xcor)
ii = length(xcor);
jj = 0;
while ii % pseudo condition
jj = jj + 1;
if(jj > ii), break; end
if(bw(ycor(jj),xcor(jj)) < 1)
xcor(jj) = [];
ycor(jj) = [];
ii = ii - 1;
end
end
end
if ~isempty(xcor)
totalElem = totalElem + numel(xcor);
for ii = 1:length(xcor)
% xcor has been calculated so far in column direction
% ycor has been calcuated in row direction
% remeber graph x direction is actually column direction
% and graph y direction is row direction
sumTotal = sumTotal + dataArray(ycor(ii),xcor(ii));
jvector_x.addElement(xcor(ii));
jvector_y.addElement(ycor(ii));
end
end
end
avg = sumTotal/totalElem;
totalElem = round(pi*a*b);
end
No comments:
Post a Comment