基于matlab的三维数据椭球拟合

Question

我正在做一个3D容积 of CT肺部图像，为了检测结节，我需要为每一个疑似结节拟合一个椭球体模型，我如何为它做一个代码？结节是怀疑是肿瘤的物体，我的算法需要检查每一个物体，并把它近似成一个椭球，从椭球参数中我们计算出8个特征来构建一个分类器，通过训练和测试数据来检测它是否是结节，所以我需要对这个椭球体进行拟合。

这里一片容积CT肺图像

​

​

在这里，另一个体积相同的切片，但它包含一个结节(黄色圆圈，有一个结节)，所以我需要我的代码来检查每个形状，以确定它是否是一个结节。

​

​

Answer 1

由于我们没有可用的3D数据集，所以我从2D开始。

首先，我们需要选择肺部，这样我们就不会计算任何其他物体，然后是里面的物体。由于这是灰度，我们首先需要以某种方式对其进行二值化。我为DIP使用了自己的DIP，这将大量使用我的growthfill，因此我强烈建议您首先阅读以下内容：

• Fracture detection in hand using image proccessing

在那里你会找到所有你需要的解释。现在为了你的任务，我会：

1. <0,765> 将RGBA转换为灰度 我只是计算强度i=R+G+B，对于24位图像，通道是8位，结果达到3*255 = 765。由于输入图像被JPEG压缩，图像中存在颜色失真和噪声，所以不要忘记这一点。
2. 裁剪出白色边框 只要投射线(扫描线)从中间的每一条边框向中间和停止，如果非白色像素击中。我用700代替765来补偿图像中的噪声。现在你得到了可用图像的边框，所以剩下的部分就会出现。
3. 计算直方图 为了补偿图像中的失真，使直方图平滑，以消除所有不必要的噪声和空白。然后从左和右找到局部最大值(red)。这将用于二进制化重持久化(中间是green) --这是我的最终结果：

​

1. 二值图像 只需将图像与直方图中的*绿色*强度相对应。因此，如果i0,i1是直方图中从左到右的局部最大强度，则对(i0+i1)/2保持不变。其结果是：

​

1. 除去除肺外的所有东西， 这很容易，只需填补黑色从外部到一些预定义的背景颜色。然后，同样的方式，所有白色的东西相邻的背景颜色。这将去除人体表面，骨骼，器官和CT机器只留下肺部。现在用预定义的肺颜色重新着色剩下的黑色。

不应该留下黑色，剩下的白色是可能的结节。

1. 处理所有剩余的白色像素 因此，只要循环通过图像和第一个白色像素命中洪水，填补它的预定义结节颜色或明确的对象索引后使用。我还区分了表面(aqua)和内部(洋红)。其结果是：

现在你可以计算每个结节的特征。如果为此编写自定义floodfill，则可以直接从它获得以下内容：

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

这一切都可以作为你的特征变量，也有助于拟合。

1. 拟合发现的曲面点 有很多方法，但我会尽可能地放松它，以提高性能和准确性。例如，您可以使用质心作为椭球中心。然后找出它的min和max距离点，并使用它们作为半轴(+/-一些正交修正)。然后在这些初始值附近搜索。有关更多信息，请参见：

- [How approximation search works](https://stackoverflow.com/a/36163847/2521214)

您将在链接的Q/A中找到使用示例。

Notes

所有的子弹都适用于3D。在构建自定义floodfill时，要小心递归尾。过多的信息会使溢出堆栈，也会大大减慢速度。在这里，我使用了很少的自定义返回参数+ growthfill来处理它，这是一个小例子：

//---------------------------------------------------------------------------
    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);
    }
//---------------------------------------------------------------------------

在这里，我使用C++代码编写了示例图像：

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;

Answer 2

您可能会对我们最近为开放源码软件Icy http://icy.bioimageanalysis.org/开发的插件感兴趣。

插件名为FitEllipsoid，它允许首先单击正交视图上的几个点，快速地将椭球与图像内容相匹配。这里有一个教程：https://www.youtube.com/watch?v=MjotgTZi6RQ

还请注意，我们在GitHub上提供了Matlab和Java源代码(但我不能提供它们，因为这是我第一次出现在网站上)。