Thread Subject: Fit sphere to planar circles

Hi all,

I'm trying to find the center and radius of a mostly spherical body in a voxel volume (femoral head in a CT scan, to be precise).

My approach so far has been to use hough circle detection in the area of interest, looped over each slice. This detection works best around the hemisphere (rather than pole) regions of the sphere, and I am using some logic to cull away any false positive circles detected by the hough transform.

My question is this: what is the best way to fit a sphere to a series of planar circles?

My data is currently in the form of:
circs = [...
  240.5275 145.0197 -564.5000 20.5078;
  241.1135 146.1915 -559.5000 22.8516;
  241.1135 146.1915 -554.5000 23.4375;
  241.6994 146.7775 -549.5000 23.4375;
  241.1135 146.7775 -544.5000 22.2656;]

The columns are the X, Y, and Z locations of 5 circle origins, and their radii.

I'm confident that the average of the X and Y columns represent the X and Y sphere centre, but I don't know the best way to get the Z location.

A simple average of Z locations is inaccurate because the detected circles are not necessarily centrally detected along the Z-direction.

Any thoughts? Perhaps I could try to fit a curve to the radii values and find its peak, but I'm not sure how to fit anything robustly to so few data points.

Any help would be very welcome.
Thanks,
Sven.

On Apr 7, 1:24=A0pm, "Sven" <sven.holco...@gmail.deleteme.com> wrote:
> Hi all,
>
> I'm trying to find the center and radius of a mostly spherical body in a =
voxel volume (femoral head in a CT scan, to be precise).
>
> My approach so far has been to use hough circle detection in the area of =
interest, looped over each slice. This detection works best around the hemi=
sphere (rather than pole) regions of the sphere, and I am using some logic =
to cull away any false positive circles detected by the hough transform.
>
> My question is this: what is the best way to fit a sphere to a series of =
planar circles?
>
> My data is currently in the form of:
> circs =3D [...
> =A0 240.5275 =A0145.0197 -564.5000 =A0 20.5078;
> =A0 241.1135 =A0146.1915 -559.5000 =A0 22.8516;
> =A0 241.1135 =A0146.1915 -554.5000 =A0 23.4375;
> =A0 241.6994 =A0146.7775 -549.5000 =A0 23.4375;
> =A0 241.1135 =A0146.7775 -544.5000 =A0 22.2656;]
>
> The columns are the X, Y, and Z locations of 5 circle origins, and their =
radii.
>
> I'm confident that the average of the X and Y columns represent the X and=
 Y sphere centre, but I don't know the best way to get the Z location.
>
> A simple average of Z locations is inaccurate because the detected circle=
s are not necessarily centrally detected along the Z-direction.
>
> Any thoughts? Perhaps I could try to fit a curve to the radii values and =
find its peak, but I'm not sure how to fit anything robustly to so few data=
 points.
>
> Any help would be very welcome.
> Thanks,
> Sven.

Sven,

I have a similar problem and I abandoned the hough transform. I assume
you have grayscale images. What I do is replicate slices to bring z
resolution ~equal to xy resolution, then I do a fft convolution with a
spherical mask (there is an fspecial3 on file exchage, you need some
normalization though) at the radius of interest. If you don't know the
radius, you can iterate over multiple radii. Take the pixel
coordinates corresponding to the max of the convolution.

Any thoughts on/need of computing the entropy of the pixels within the
sphere without requiring a supercomputing center? I'm actively
researching this. Let me know if you want to chat more about this.

Thanks,
Trance.

trancemissionxxi@gmail.com wrote in message <c91073f0-2621-4779-b9a0-b2d3ad8083d2@z1g2000yqn.googlegroups.com>...
> On Apr 7, 1:24=A0pm, "Sven" <sven.holco...@gmail.deleteme.com> wrote:
> > Hi all,
> >
> > I'm trying to find the center and radius of a mostly spherical body in a =
> voxel volume (femoral head in a CT scan, to be precise).
> >
> > My approach so far has been to use hough circle detection in the area of =
> interest, looped over each slice. This detection works best around the hemi=
> sphere (rather than pole) regions of the sphere, and I am using some logic =
> to cull away any false positive circles detected by the hough transform.
> >
> > My question is this: what is the best way to fit a sphere to a series of =
> planar circles?
> >
> > My data is currently in the form of:
> > circs =3D [...
> > =A0 240.5275 =A0145.0197 -564.5000 =A0 20.5078;
> > =A0 241.1135 =A0146.1915 -559.5000 =A0 22.8516;
> > =A0 241.1135 =A0146.1915 -554.5000 =A0 23.4375;
> > =A0 241.6994 =A0146.7775 -549.5000 =A0 23.4375;
> > =A0 241.1135 =A0146.7775 -544.5000 =A0 22.2656;]
> >
> > The columns are the X, Y, and Z locations of 5 circle origins, and their =
> radii.
> >
> > I'm confident that the average of the X and Y columns represent the X and=
> Y sphere centre, but I don't know the best way to get the Z location.
> >
> > A simple average of Z locations is inaccurate because the detected circle=
> s are not necessarily centrally detected along the Z-direction.
> >
> > Any thoughts? Perhaps I could try to fit a curve to the radii values and =
> find its peak, but I'm not sure how to fit anything robustly to so few data=
> points.
> >
> > Any help would be very welcome.
> > Thanks,
> > Sven.
>
> Sven,
>
> I have a similar problem and I abandoned the hough transform. I assume
> you have grayscale images. What I do is replicate slices to bring z
> resolution ~equal to xy resolution, then I do a fft convolution with a
> spherical mask (there is an fspecial3 on file exchage, you need some
> normalization though) at the radius of interest. If you don't know the
> radius, you can iterate over multiple radii. Take the pixel
> coordinates corresponding to the max of the convolution.
>
> Any thoughts on/need of computing the entropy of the pixels within the
> sphere without requiring a supercomputing center? I'm actively
> researching this. Let me know if you want to chat more about this.
>
> Thanks,
> Trance.

Thanks for the input, Trance.
I can see how the 3d convolution would simplify the hough transform approach to essentially a single operation. I wasn't able to find an fspecial3 entry in the file exchange via any expected keywords... do you have a direct link?

I don't have any experience with computing pixel entropy, so I'm afraid I'm not going to be of much assistance there yet.

Without having yet implemented your spherical fft solution, I'm wondering how well it has worked for you. I'm not sure if you're doing the exact same thing (femoral head from CT), but if you are, perhaps you have encountered the problem of the acetabular surrounding the femoral head being basically spherical too, just with a slightly larger radius.

I'm thinking that it might even be a good idea to go once through with a slice by slice circle detection from hough transforms, then produce a new voxel volume with pixels on the boundaries of all detected circles turned on, and try your 3d fft on *that* volume, rather than the original grayscale volume.

Cheers,
Sven.

Sven,

For some reason fspecial3 is not showing up on file exchange anymore.
I googled it and the link is:

http://www.mathworks.com/matlabcentral/files/15981/fspecial3.m

I do 3D cell nuclei tracking, which are almost spherical, but I know
about the acetabular issue from papers I read on hough transform
(tried using it before for my problem).

Indeed the fft can become computationally intensive if the image and
mask are large (for me largest volume element is 49x49x49 pixels).
I like your idea on doing the convolution on the hough transformed
circles/disks. Since you are done with using the grayscale and have
binary data for the slice circles, why not resize the image and do
convolution on smaller image size (if large matrices are a problem).
Alternatively, you can use a scanning box algorithm on your new voxel
volume with detected disks.

There's also this paper (not sure if you have access):
http://dx.doi.org/10.1016/j.patrec.2005.11.019

Cheers.

trancemissionxxi@gmail.com wrote in message <11b5e1ed-fdcf-4eea-a9d5-f7c49b221c40@a7g2000yqk.googlegroups.com>...
> Sven,
>
> For some reason fspecial3 is not showing up on file exchange anymore.
> I googled it and the link is:
>
> http://www.mathworks.com/matlabcentral/files/15981/fspecial3.m
>
> I do 3D cell nuclei tracking, which are almost spherical, but I know
> about the acetabular issue from papers I read on hough transform
> (tried using it before for my problem).
>
> Indeed the fft can become computationally intensive if the image and
> mask are large (for me largest volume element is 49x49x49 pixels).
> I like your idea on doing the convolution on the hough transformed
> circles/disks. Since you are done with using the grayscale and have
> binary data for the slice circles, why not resize the image and do
> convolution on smaller image size (if large matrices are a problem).
> Alternatively, you can use a scanning box algorithm on your new voxel
> volume with detected disks.
>
> There's also this paper (not sure if you have access):
> http://dx.doi.org/10.1016/j.patrec.2005.11.019
>
> Cheers.

Thanks Trance, that paper actually helped quite a lot. It essentially has an implementation of my "find the best sphere to fit a set of planar circles" problem, and it ends up boiling down to pretty simple trigonometry after you make the (fairly reasonable) assumption that the circles all have the same (x,y) location.

I had a little attempt at the fspecial3 approach, but as you know it is a little computationally expensive, and I'm quite happy with the results so far from the hough circles approach. The CT scans are often higher than 49x49 pixels across the area of interest without downsampling, so hough circles were much faster to compute than the full fft.

One further question to do with the next step of what I'm considering. Do you know or use any kind of radial filtering? I'm wondering if I could somehow filter the image volume after I find the femoral head using the head centre and radius, in a way that would highlight the acetabular surface. I'm pretty sure a simple threshold outside the sphere would do a reasonable job, but I'd just like to know if what I'm thinking is effective and maybe even commonly used.

Cheers,
Sven.

"Sven" <sven.holcombe@gmail.deleteme.com> wrote in message <gro7to$k0a$1@fred.mathworks.com>...
> trancemissionxxi@gmail.com wrote in message <11b5e1ed-fdcf-4eea-a9d5-f7c49b221c40@a7g2000yqk.googlegroups.com>...
> > Sven,
> >
> > For some reason fspecial3 is not showing up on file exchange anymore.
> > I googled it and the link is:
> >
> > http://www.mathworks.com/matlabcentral/files/15981/fspecial3.m
> >
> > I do 3D cell nuclei tracking, which are almost spherical, but I know
> > about the acetabular issue from papers I read on hough transform
> > (tried using it before for my problem).
> >
> > Indeed the fft can become computationally intensive if the image and
> > mask are large (for me largest volume element is 49x49x49 pixels).
> > I like your idea on doing the convolution on the hough transformed
> > circles/disks. Since you are done with using the grayscale and have
> > binary data for the slice circles, why not resize the image and do
> > convolution on smaller image size (if large matrices are a problem).
> > Alternatively, you can use a scanning box algorithm on your new voxel
> > volume with detected disks.
> >
> > There's also this paper (not sure if you have access):
> > http://dx.doi.org/10.1016/j.patrec.2005.11.019
> >
> > Cheers.
>
> Thanks Trance, that paper actually helped quite a lot. It essentially has an implementation of my "find the best sphere to fit a set of planar circles" problem, and it ends up boiling down to pretty simple trigonometry after you make the (fairly reasonable) assumption that the circles all have the same (x,y) location.
>
> I had a little attempt at the fspecial3 approach, but as you know it is a little computationally expensive, and I'm quite happy with the results so far from the hough circles approach. The CT scans are often higher than 49x49 pixels across the area of interest without downsampling, so hough circles were much faster to compute than the full fft.
>
> One further question to do with the next step of what I'm considering. Do you know or use any kind of radial filtering? I'm wondering if I could somehow filter the image volume after I find the femoral head using the head centre and radius, in a way that would highlight the acetabular surface. I'm pretty sure a simple threshold outside the sphere would do a reasonable job, but I'd just like to know if what I'm thinking is effective and maybe even commonly used.
>

Just a little update if anyone's interested.
I found that the above radial filter worked quite well in my application. Let me explain:

The hough circle detection algorithm I was using was "houghcircles" from the file exchange. By default, this file uses a 'canny' filter on a given image. A minor problem I was having was that the femoral head and acetabular are in close proximity, and it was difficult to distinguish the edge of one from the edge of the other.
So, I went for a 2 phase approach.

First phase:
- Simple circle detection (as per houghcircles) in order to find a *pretty good* [x0,y0] location of the centre of my sphere. (the paper linked above has some tips)

Second phase:
- A second round of hough circle detection, however this time I don't use the raw input image. Instead, something like the following:
% xVoxels,yVoxels : vectors for x and y pixel locations
% im : raw image at a given height (size matches yVoxels/xVoxels)

[xPx, yPx] = meshgrid(xVoxels, yVoxels); % Cartesian coords of pixel locs
[th,rho] = cart2pol(xPx-x0, yPx-y0);
[px,py] = gradient(im,xVoxels,yVoxels); % XY Gradient filter
pr = sin(th).*py + cos(th).*px; % Radial gradient filter
bw_filt = bwmorph(pr<-50, 'skel')

This gives a thin line at the boundaries where the gradient of the image is negative 50 or lower (in my case, this fits well to the gradient of bone to non-bone).

I then provide this bw_filt image instead of the raw image to the houghcircles function (and I modified the function so it didn't use a canny filter on this image). The result was a much improved detection of the outer femoral head surface, since it effectively went from detecting:
femoral head inner edge
femoral head outer edge
acetabular outer edge
acetabular inner edge

to just femoral head outer edge and acetabular inner edge.

Anyway, just thought that might come in handy.

Cheers,
Sven.