I have what amounts to an xray image of a circular xray field - the xray beam is a cone. The beam is supposed to be aimed at a metal ball with the shadow of the ball hopefully centered on the xray beam. So you basically have two concentric circles with inner one a lighter shade of grey and the outer basically black (the exposed part of the film). I hope you're still with me.

What I want to do is (1) find the centers of both circles and (2) find the eccentricity of the two circles. My first approach was to try to find the boundary of each circle and then perform some measurements on those. However, today I began thinking maybe I could simply locate two points on each circle, form a chord, and determine the center. This assumes a perfect circle, which may not be the case. I could do this several times and calculate a standard deviation.

Now to my question. Which way would you pros attack this. Perhaps there is a completely different way I need to look at this.

Thanks, Jeff

doc regionprops

"Jeff" wrote in message <ilvdep$1tg$1@ginger.mathworks.com>...
> I have what amounts to an xray image of a circular xray field - the xray beam is a cone. The beam is supposed to be aimed at a metal ball with the shadow of the ball hopefully centered on the xray beam. So you basically have two concentric circles with inner one a lighter shade of grey and the outer basically black (the exposed part of the film). I hope you're still with me.
> 
> What I want to do is (1) find the centers of both circles and (2) find the eccentricity of the two circles. My first approach was to try to find the boundary of each circle and then perform some measurements on those. However, today I began thinking maybe I could simply locate two points on each circle, form a chord, and determine the center. This assumes a perfect circle, which may not be the case. I could do this several times and calculate a standard deviation.
> 
> Now to my question. Which way would you pros attack this. Perhaps there is a completely different way I need to look at this.
> 
> Thanks, Jeff

To specify a circle we need at least three points. So, you find three points corresponding to the circles and then find the respective centres. Then, you can easily find the radius of the circles. The one having the larger radius will be the outer circle and the one with the smaller radius will be the inner circle. If you do not know know how to compute center of a circle from three points on its circumference, you look search the term "Center of Points" + Yumnam, that will point to a website I wrong long time back describing it, but I could not remember the site. 

Yumnam Kirani Singh
Tronglaobi Awang Leikai

"Jeff" wrote in message <ilvdep$1tg$1@ginger.mathworks.com>...
> I have what amounts to an xray image of a circular xray field - the xray beam is a cone. The beam is supposed to be aimed at a metal ball with the shadow of the ball hopefully centered on the xray beam. So you basically have two concentric circles with inner one a lighter shade of grey and the outer basically black (the exposed part of the film). I hope you're still with me.
> 
> What I want to do is (1) find the centers of both circles and (2) find the eccentricity of the two circles. My first approach was to try to find the boundary of each circle and then perform some measurements on those. However, today I began thinking maybe I could simply locate two points on each circle, form a chord, and determine the center. This assumes a perfect circle, which may not be the case. I could do this several times and calculate a standard deviation.
> 
> Now to my question. Which way would you pros attack this. Perhaps there is a completely different way I need to look at this.
> 
> Thanks, Jeff
Happy Holi!
You can visit the site http://www.oocities.org/kiranisingh/center.html and use formula 2 to find the center of circles from three points.  The three points should corresponds to the circumference of the circle for which you want to compute center.  

Yumnam Kirani Singh
Tronglaobi Awang Leikai

"Yumnam Kirani" <kirani.singh@gmail.com> wrote in message <ilvilk$epp$1@ginger.mathworks.com>...
>
> 
> To specify a circle we need at least three points. So, you find three points corresponding to the circles and then find the respective centres. 
==============

Not a good idea. There will always be errors in identifying the location of points on the boundary of a shape. Therefore using only 3 data points will make the result sensitivie to errors. REGIONPROPS uses all boundary points to derive a much more exact result.

REGIONPROPS requires the Image Processing Toolbox. If you don't have that, there is my ellipsefit function below. A number of alternative shape fitting functions are also available on the File Exchange.



function report=ellipsefit(XY)
%ELLIPSEFIT - form 2D ellipse fit to given x,y data
%
%  report=ellipsefit(XY)
%
%in:
%
% XY: Input matrix of 2D coordinates to be fit. Each column XY(:,i) is [xi;yi]
%
%out: Finds the ellipse fitting the input data parametrized both as
%     A*x^2+B*x*y C*y^2+D*x+E*y=1 and [x-x0,y-y0]*Q*[x-x0;y-y0]=1
%
% report: a structure output with the following fields
%    
%    report.Q: the matrix Q
%    report.d: the vector [x0,y0]
%    report.ABCDE: the vector [A,B,C,D,E]
%    report.AxesDiams: The minor and major ellipse diameters
%    report.theta: The counter-clockwise rotation of the ellipse.
%
%NOTE: The code will give errors if the data fit traces out a non-elliptic or
%      degenerate conic section.
%
%See also ellipsoidfit

X=XY(1,:).';
Y=XY(2,:).';

M= [X.^2, X.*Y, Y.^2, X, Y, -ones(size(X,1),1)];

[U,S,V]=svd(M,0);
ABCDEF=V(:,end);

if size(ABCDEF,2)>1
    
    error 'Data cannot be fit with unique ellipse'
else

    ABCDEF=num2cell(ABCDEF);
end

[A,B,C,D,E,F]=deal(ABCDEF{:});


Q=[A, B/2;B/2 C];
x0=-Q\[D;E]/2;
dd=F+x0'*Q*x0;

Q=Q/dd;

[R,eigen]=eig(Q);
eigen=eigen([1,4]);

if ~all(eigen>=0), error 'Fit produced a non-elliptic conic section'; end

idx=eigen>0;
eigen(idx)=1./eigen(idx);
AxesDiams = 2*sqrt(eigen);

theta=atand(tand(-atan2(R(1),R(2))*180/pi));


report.Q=Q;
report.d=x0(:).';
report.ABCDE=[A, B, C, D, E]/F;
report.AxesDiams=sort(AxesDiams(:)).';
report.theta=theta;

"Matt J" wrote in message <ilvgno$9tp$1@ginger.mathworks.com>...
> 
> doc regionprops

Thanks, Matt. I'll take a look at that.

Thanks, Yumnan. Although I haven't looked at the formlae you provided, the theory I was using is that two points on the circle form a chord. A perpendicular bisector of the chord will yield a diameter, the center of which is the center of the circle.

Matt, you are correct in that this method can yield variable reuslts. My initial thought was to repeat this process several times providing a number of solutions. Hopefully the average would yield an accurate result - a porr man's Monte carlo if you will.

"Jeff" wrote in message <ilvka5$jrp$1@ginger.mathworks.com>...
>
> Matt, you are correct in that this method can yield variable reuslts. My initial thought was to repeat this process several times providing a number of solutions. Hopefully the average would yield an accurate result - a porr man's Monte carlo if you will.
========================

It wouldn't be bad. But if the center is the only thing you want, the optimal poor man's Monte Carlo method would be to compute the centroid directly. Here is an example for a square, but it would work for any uniform shape.

>> u=[0 1 1 1 0 0 0 ]; v=[ 0 0 1 1 1 0 0]; Image=u'*v

Image =

     0     0     0     0     0     0     0
     0     0     1     1     1     0     0
     0     0     1     1     1     0     0
     0     0     1     1     1     0     0
     0     0     0     0     0     0     0
     0     0     0     0     0     0     0
     0     0     0     0     0     0     0


>> [x,y]=ndgrid(1:length(Image));

>> Xcenter=mean(x(Image>0))

Xcenter =

     3

>> Ycenter=mean(y(Image>0))

Ycenter =

     4

Jeff:
Is the region outside the outer dark circle also light?  So basically
we have light outer stuff (from dark ring to edges of image that
didn't get any X rays, then the dark ring which had max X rays, then
the light center disc which had lesser x-ray exposure).  Is that
correct?  It would have been helpful if you had uploaded an image.

To find the center you'd basically threshold, label, and then use
regionprops to find the weighted centroid and the eccentricity.  I do
that same thing in my BlobsDemo tutorial, where I find the centroids
of light coins on a dark background.

http://www.mathworks.com/matlabcentral/fileexchange/25157
It's all there in well commented, excruciating detail that any novice
should be able to follow.

You'd do the same thing except that I'd recommend you use
"WeightedCentroid" instead of "Centroid" for the light center ball ,
but you'd use "Centroid" for the outer dark ring.

For the outer dark circle, you'd have to threshold to get the dark
ring, instead of the bright circle.  I don't think you'd need to do a
hole fill to get the centroid.  Even though the centroid is outside of
the dark ring, it should still get the centroid in the center of the
dark ring.

That's how I'd attack it, and I guess you'd call me a pro - everyone
else does and I've been doing this over 30 years, and got my Ph.D. in
radiological image processing and still do some radiological imaging.

-ImageAnalyst

ImageAnalyst,

You are correct. The image has a light outer area with a dark annulus and light inner circle. As you point out the dark middle annulus is a result of maximum x-ray exposure. The piece of film is square. I'll try to find a free hosting site so I can post a pic. Can you recommend a good free site to put my image in?

Interesting that you were involved in radiological imaging. I'm a medical physicist in a radiation therapy center. The project I'm working on is for quality assurance tests. 

Jeff

Here is a link to a sample of the types of images  I am trying to process. Keep in mind, the density and annotations may change from day to day. So again to recap, the dark outer ring is a circular x-ray beam aimed at a metal ball which is demonstrated by the light inner area. The assumption is that the x-ray beam and ball are coincident - that is what we are checking. Also the film is perpendicular to the x-ray beam. If it isn't then the beam will not be a perfect circle and will have an eccentricity not equal to 1.0.

Link: http://img196.imageshack.us/img196/1627/aqabase001.png

Thanks again for everyone's help.
Jeff

My original tif image is loaded as a uint16 data type. Do I need to convert this to something else?

Jeff

On Mar 19, 5:33=A0pm, "Jeff " <physic...@bellsouth.net> wrote:
> My original tif image is loaded as a uint16 data type. Do I need to conve=
rt this to something else?
>
> Jeff
---------------------------------------------------------------------------=
--------
Jeff:
I can't download that image.  Right clicking only allows me to save it
as an htm file, which is no good.  Try uploadhouse.com.

You don't need to convert it to something else.  MATLAB works with 16
bit images.
I'll see if I can modify me BlobsDemo script later if I have time (and
the image.)
ImageAnalyst

Jeff:
Try this code:

% Demo find the weighted centroid and eccentricity of a circle and a
ring.
% By ImageAnalyst
%
% IMPORTANT: The newsreader may break long lines into multiple lines.
% Be sure to join any long lines that got split into multiple single
lines.
% These can be found by the red lines on the left side of your
% text editor, which indicate syntax errors, or else just run the
% code and it will stop at the split lines with an error.
% function test()
% Change the current folder to the folder of this m-file.
if(~isdeployed)
        cd(fileparts(which(mfilename)));
end
clc;    % Clear command window.
clear;  % Delete all variables.
close all;      % Close all figure windows except those created by
imtool.
imtool close all;       % Close all figure windows created by imtool.
workspace;      % Make sure the workspace panel is showing.
fontSize = 16;
% Read in a standard MATLAB color demo image.
folder = 'C:\Documents and Settings\user\My Documents\Downloads';
baseFileName = 'xray_ball.png';
% Get the full filename, with path prepended.
fullFileName = fullfile(folder, baseFileName);
if ~exist(fullFileName, 'file')
	% Didn't find it there.  Check the search path for it.
	fullFileName = baseFileName; % No path this time.
	if ~exist(fullFileName, 'file')
		% Still didn't find it.  Alert user.
		errorMessage = sprintf('Error: %s does not exist.', fullFileName);
		uiwait(warndlg(errorMessage));
		return;
	end
end
grayImage = imread(fullFileName);
% Get the dimensions of the image.  numberOfColorBands should be = 1.
[rows columns numberOfColorBands] = size(grayImage);
if numberOfColorBands > 1
	grayImage = rgb2gray(grayImage);
end
% Display the original gray scale image.
subplot(2, 3, 1);
imshow(grayImage, []);
axis on;
title('Original Grayscale Image', 'FontSize', fontSize);
% Enlarge figure to full screen.
set(gcf, 'Position', get(0,'Screensize'));
set(gcf,'name','Demo by ImageAnalyst','numbertitle','off')

% imtool(grayImage);

% Let's compute and display the histogram.
[pixelCount grayLevels] = imhist(grayImage);
subplot(2, 3, 2);
bar(pixelCount);
title('Histogram of original image', 'FontSize', fontSize);
xlim([0 grayLevels(end)]); % Scale x axis manually.

% Find the central ball
thresholdValue = 67;
% Threshold the image to get a binary image (only 0's and 1's) of
class "logical."
binaryImage = grayImage > thresholdValue; % Bright objects will be the
chosen if you use >.
% Get rid of bright outer stuff.
binaryImage = imclearborder(binaryImage);
% Do a "hole fill" to get rid of any background pixels inside the
blobs.
binaryImage = imfill(binaryImage, 'holes');
% Get rid of small things.
binaryImage = bwareaopen(binaryImage, 1000);
subplot(2, 3, 3);
imshow(binaryImage, []);
axis on;
title('Binary Image of the Ball', 'FontSize', fontSize);
% Label the image and make the measurements
labeledImage = bwlabel(binaryImage, 8);     % Label each blob so we
can make measurements of it
% Get all the blob properties.  Can only pass in grayImage in version
R2008a and later.
blobMeasurements = regionprops(labeledImage, grayImage,
'WeightedCentroid', 'Eccentricity');
ballCentroid = blobMeasurements(1).WeightedCentroid;		% Get centroid.
ballEccentricity = blobMeasurements(1).Eccentricity;
fprintf(1, 'The weighted centroid is at (%.3f, %.3f)\n',
ballCentroid(1), ballCentroid(2));
fprintf(1, 'The eccentricity of the ball = %.5f\n', ballEccentricity);
% Find the ring.
thresholdValue = 66;
% Threshold the image to get a binary image (only 0's and 1's) of
class "logical."
binaryRing = grayImage < thresholdValue; % Bright objects will be the
chosen if you use >.
% Get rid of bright outer stuff.
binaryRing = imclearborder(binaryRing);
% Do a "hole fill" to get rid of any background pixels inside the
blobs.
binaryRing = imfill(binaryRing, 'holes');
% Get rid of small things.
binaryRing = bwareaopen(binaryRing, 1000);
subplot(2, 3, 4);
imshow(binaryRing, []);
axis on;
title('Binary Image of the Ring', 'FontSize', fontSize);
% Label the image and make the measurements
labeledImage = bwlabel(binaryRing, 8);     % Label each blob so we can
make measurements of it
% Get all the blob properties.  Can only pass in grayImage in version
R2008a and later.
blobMeasurements = regionprops(labeledImage, grayImage,
'WeightedCentroid', 'Eccentricity');

ringCentroid = blobMeasurements(1).WeightedCentroid;		% Get centroid.
ringEccentricity = blobMeasurements(1).Eccentricity;
fprintf(1, 'The weighted centroid is at (%.3f, %.3f)\n',
ringCentroid(1), ringCentroid(2));
fprintf(1, 'The eccentricity of the ring = %.5f\n', ringEccentricity);

msgbox('Done with demo.  Look in the command window');

IA,

Thanks for the code. Can't wait to work on it tomorrow. I especially liked the code to remove the outer bright sections and the small areas (noise). Guess I could have read through more of the image processing toolbox user guide. 

One question I have after looking at your code. How did you choose a threshold value of 66 and 67 when selecting the inner/outer ball without knowing what the actual grey values are?when I have run my code on different films of the same setup I have had to use different threshold values. 

Thanks,
Jeff

Jeff:
That's correct.  I just used imtool to determine what threshold values
to use.  I just moved the cursor over the edge of the ball and ring
and found the gray level that was at what I thought was the separation
circle, so to speak.
ImageAnalyst