Ellipsoid fitting for 3D data on matlab

Question

I am working on a 3D volume of CT lung images, In order to detect nodules I need to fit an ellipsoid model for each suspected nodule, How can I make a code for that ??? Nodule is the suspected object to be a tumor, my algorithm needs to check every object, and approximate it to an ellipsoid, and from the ellipsoid parameters we calculate 8 features to build a classifier which detects whether it a nodule or not through training and testing data, so I need to fit such ellipsoid

here one slice of the volume CT lung image

here another slice of the same volume but it contains a nodule (the yellow circle there is a nodule) so I need my code to check every shape to determine whether it is a nodule or not

Answer 1

you may be interested in a recent plugin that we developed for the open-source software Icy http://icy.bioimageanalysis.org/

The plugin name is FitEllipsoid, it allows rapidly fitting an ellipsoid to the image contents by first clicking on a few points on the orthogonal view. A tutorial is available here: https://www.youtube.com/watch?v=MjotgTZi6RQ

Also note that we provide Matlab and Java source codes on GitHub (but I cannot provide them since it is my first appearance on the website).

Answer 2

As we do not have 3D dataset at disposal I start with 2D.

So first we need to select the lungs so we do not count any other objects then those inside. As this is gray scale we first need to binarize it somehow. I use my own class picture for DIP and this will heavily use my growthfill so I strongly recommend to first read this:

• Fracture detection in hand using image proccessing

Where you will find all explanations you need for this. Now for your task I would:

1. turn RGBA to grayscale <0,765>

   I just compute intensity i=R+G+B as for 24 bit image the channels are 8 bit the result is up to 3*255 = 765. As the input image was compressed by JPEG there are distortions in color and also noise in the image so do not forget about that.

2. crop out the white border

   Just cast rays (scan lines) from middle of each outer line of the border towards middle and stop if non white-ish pixel hit. I treshold with 700 instead of 765 to compensate the noise in image. Now you got the bounding box of usable image so crop out the rest.

3. compute histogram

   To compensate distortions in image smooth the histogram enough to remove all unwanted noise and gaps. Then find the local maximum from left and from right (red). This will be used for binarisation treshold (middle between them green) This is my final result:

   

4. binarise image

   Just treshold image against the **green* intensity from histogram. So if i0,i1 are the local maximum intensities from left and right in the histogram then treshold against (i0+i1)/2. This is the result:

   

5. remove everything except lungs

   That is easy just fill the black from outside to some predefined background color. Then the same way all the white stuff neighboring background color. That will remove human surface, skeleton, organs and the CT machine leaving just the lungs. Now recolor the remaining black with some predefined Lungs color.

   

   There should be no black color left and the remaining white are the possible nodules.

6. process all remaining white pixels

   So just loop through the image and on first white pixel hit flood fill it with predefined nodule color or distinct object index for latter use. I also distinct the surface (aqua) and the inside (magenta). This is the result:

   

   Now you can compute your features per nodule. If you code your custom floodfill for this then you can obtain from it directly things like:

   • Volume in [pixels]
   • Surface area in [pixels]
   • Bounding box
   • Position (relative to Lungs)
   • Center of gravity or centroid

   Which all can be used as your feature variables and also to help with fitting.

7. fit the found surface points

   There are many methods for this but I would ease up it as much as I could to improve performance and accuracy. For example you can use centroid as your ellipsoid center. Then find the min and max distant points from it and use them as semi-axises (+/- some orthogonality corrections). And then just search around these initial values. For more info see:

   • How approximation search works

   You will find examples of use in linked Q/As there.

[Notes]

All of the bullets are applicable in 3D. While constructing custom floodfill be careful with the recursion tail. Too much info will really fast overflow your stack and also slow things down considerably. Here small example of how I deal with it with few custom return parameters + growthfill I used:

//---------------------------------------------------------------------------
    void growfill(DWORD c0,DWORD c1,DWORD c2);                      // grow/flood fill c0 neigbouring c1 with c2
    void floodfill(int x,int y,DWORD c);                            // flood fill from (x,y) with color c
    DWORD _floodfill_c0,_floodfill_c1;                              // recursion filled color and fill color
    int   _floodfill_x0,_floodfill_x1,_floodfill_n;                 // recursion bounding box and filled pixel count
    int   _floodfill_y0,_floodfill_y1;
    void  _floodfill(int x,int y);                                  // recursion for floodfill
//---------------------------------------------------------------------------
void picture::growfill(DWORD c0,DWORD c1,DWORD c2)
    {
    int x,y,e;
    for (e=1;e;)
     for (e=0,y=1;y<ys-1;y++)
      for (   x=1;x<xs-1;x++)
       if (p[y][x].dd==c0)
        if ((p[y-1][x].dd==c1)
          ||(p[y+1][x].dd==c1)
          ||(p[y][x-1].dd==c1)
          ||(p[y][x+1].dd==c1)) { e=1; p[y][x].dd=c2; }
    }
//---------------------------------------------------------------------------
void picture::_floodfill(int x,int y)
    {
    if (p[y][x].dd!=_floodfill_c0) return;
    p[y][x].dd=_floodfill_c1;
    _floodfill_n++;
    if (_floodfill_x0>x) _floodfill_x0=x;
    if (_floodfill_y0>y) _floodfill_y0=y;
    if (_floodfill_x1<x) _floodfill_x1=x;
    if (_floodfill_y1<y) _floodfill_y1=y;
    if (x>   0) _floodfill(x-1,y);
    if (x<xs-1) _floodfill(x+1,y);
    if (y>   0) _floodfill(x,y-1);
    if (y<ys-1) _floodfill(x,y+1);
    }
void picture::floodfill(int x,int y,DWORD c)
    {
    if ((x<0)||(x>=xs)||(y<0)||(y>=ys)) return;
    _floodfill_c0=p[y][x].dd;
    _floodfill_c1=c;
    _floodfill_n=0;
    _floodfill_x0=x;
    _floodfill_y0=y;
    _floodfill_x1=x;
    _floodfill_y1=y;

    _floodfill(x,y);
    }
//---------------------------------------------------------------------------

And here C++ code I did the example images with:

picture pic0,pic1;
    // pic0 - source img
    // pic1 - output img
int x0,y0,x1,y1,x,y,i,hist[766];
color c;
// copy source image to output
pic1=pic0;
pic1.pixel_format(_pf_u);           // grayscale <0,765>
//                    0xAARRGGBB
const DWORD col_backg=0x00202020;   // gray
const DWORD col_lungs=0x00000040;   // blue
const DWORD col_out  =0x0000FFFF;   // aqua nodule surface
const DWORD col_in   =0x00800080;   // magenta nodule inside
const DWORD col_test =0x00008040;   // green-ish distinct color just for safe recoloring

// [remove white background]
// find white background area (by casting rays from middle of border into center of image till non white pixel hit)
for (x0=0        ,y=pic1.ys>>1;x0<pic1.xs;x0++) if (pic1.p[y][x0].dd<700) break;
for (x1=pic1.xs-1,y=pic1.ys>>1;x1>      0;x1--) if (pic1.p[y][x1].dd<700) break;
for (y0=0        ,x=pic1.xs>>1;y0<pic1.ys;y0++) if (pic1.p[y0][x].dd<700) break;
for (y1=pic1.ys-1,x=pic1.xs>>1;y1>      0;y1--) if (pic1.p[y1][x].dd<700) break;
// crop it away
pic1.bmp->Canvas->Draw(-x0,-y0,pic1.bmp);
pic1.resize(x1-x0+1,y1-y0+1);

// [prepare data]
// raw histogram
for (i=0;i<766;i++) hist[i]=0;
for (y=0;y<pic1.ys;y++)
 for (x=0;x<pic1.xs;x++)            
  hist[pic1.p[y][x].dd]++;
// smooth histogram a lot (remove noise and fill gaps due to compression and scanning nature of the image)
for (x=0;x<100;x++)
    {
    for (i=0;i<765;i++) hist[i]=(hist[i]+hist[i+1])>>1;
    for (i=765;i>0;i--) hist[i]=(hist[i]+hist[i-1])>>1;
    }
// find peaks in histogram (for tresholding)
for (x=0,x0=x,y0=hist[x];x<766;x++)
    {
    y=hist[x];
    if (y0<y) { x0=x; y0=y; }
    if (y<y0>>1) break;
    }
for (x=765,x1=x,y1=hist[x];x>=0;x--)
    {
    y=hist[x];
    if (y1<y) { x1=x; y1=y; }
    if (y<y1>>1) break;
    }
// binarize image (tresholding)
i=(x0+x1)>>1;                       // treshold with middle intensity between peeks
pic1.pf=_pf_rgba;                   // result will be RGBA
for (y=0;y<pic1.ys;y++)
 for (x=0;x<pic1.xs;x++)
  if (pic1.p[y][x].dd>=i) pic1.p[y][x].dd=0x00FFFFFF;
   else                   pic1.p[y][x].dd=0x00000000;
pic1.save("out0.png");
// recolor outer background
for (x=0;x<pic1.xs;x++) pic1.p[        0][x].dd=col_backg;  // render rectangle along outer border so the filling starts from there
for (x=0;x<pic1.xs;x++) pic1.p[pic1.ys-1][x].dd=col_backg;
for (y=0;y<pic1.ys;y++) pic1.p[y][        0].dd=col_backg;
for (y=0;y<pic1.ys;y++) pic1.p[y][pic1.xs-1].dd=col_backg;
pic1.growfill(0x00000000,col_backg,col_backg);              // fill its black content outside in
// recolor white human surface and CT machine
pic1.growfill(0x00FFFFFF,col_backg,col_backg);
// recolor Lungs dark matter
pic1.growfill(0x00000000,col_backg,col_lungs);              // outer border
pic1.growfill(0x00000000,col_lungs,col_lungs);              // fill its black content outside in
pic1.save("out1.png");
// find/recolor individual nodules
for (y=0;y<pic1.ys;y++)
 for (x=0;x<pic1.xs;x++)
  if (pic1.p[y][x].dd==0x00FFFFFF)
    {
    pic1.floodfill(x,y,col_test);
    pic1.growfill(col_lungs,col_test,col_out);
    pic1.floodfill(x,y,col_in);
    }
pic1.save("out2.png");

// render histogram
for (x=0;(x<766)&&(x>>1<pic1.xs);x++)
 for (y=0;(y<=hist[x]>>6)&&(y<pic1.ys);y++)
   pic1.p[pic1.ys-1-y][x>>1].dd=0x000040FF;
for (x=x0        ,y=0;(y<=100)&&(y<pic1.ys);y++) pic1.p[pic1.ys-1-y][x>>1].dd=0x00FF0000;
for (x=x1        ,y=0;(y<=100)&&(y<pic1.ys);y++) pic1.p[pic1.ys-1-y][x>>1].dd=0x00FF0000;
for (x=(x0+x1)>>1,y=0;(y<=100)&&(y<pic1.ys);y++) pic1.p[pic1.ys-1-y][x>>1].dd=0x0000FF00;