matlab - 在Matlab上用于3D数据的椭圆拟合

Question

我正在研究 CT 肺图像的 3D体积，为了检测结节，我需要为每个可疑结节拟合一个椭圆体模型，我该如何为此编写代码？
结节是怀疑为肿瘤的物体，我的算法需要检查每个物体，并将其近似为椭球，然后根据椭球参数计算8个特征，以建立分类器，该分类器通过训练和测试来检测是否为结节数据，所以我需要适合这样的椭球

这是一张容积的 CT肺图像

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

Answer 1

由于我们没有 3D 数据集，因此我从 2D 开始。

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

Fracture detection in hand using image proccessing

您将在这里找到所需的所有说明。现在执行您的任务:

将RGBA转换为灰度<0,765>

我只是针对24位图像计算强度i=R+G+B，通道为8位，结果取决于3*255 = 765。由于输入图像是通过JPEG压缩的，因此图像中会出现颜色失真和噪点，因此请不要忘记这一点。

裁剪出白色边框

只需将光线(扫描线)从边框的每个外线的中间投射到中间，并在非白色像素命中时停止。我使用700而不是765来补偿图像中的噪点。现在，您获得了可用图像的边界框，因此将其余部分裁剪掉。

计算直方图

为了补偿图像中的失真，请对直方图进行足够的平滑处理，以消除所有不必要的噪声和间隙。然后从左侧和右侧找到本地最大值(红色)。这将用于二值化阈值(中间是绿色)，这是我的最终结果:

二进制图像

只需将图像与直方图的“绿色”强度相对应即可。因此，如果i0,i1是直方图中左侧和右侧的局部最大强度，则针对(i0+i1)/2设置阈值。结果如下:

删除除肺部以外的所有内容

只需从外部将黑色填充为某些预定义的背景颜色即可。然后以相同的方式将所有白色填充物邻近背景色。这将去除人体表面，骨骼，器官和 CT 机器，仅留下肺部。现在，使用一些预定义的“肺”颜色为其余的黑色重新着色。



不应留有黑色，剩余的白色是可能的结节。

处理所有剩余的白色像素

因此，只需循环遍历图像，并在第一个白色像素上命中洪水，即可为其填充预定义的结节颜色或不同的对象索引，以备后用。我还区分了表面( aqua )和内部(洋红色)。结果如下:



现在，您可以计算每个结节的特征。如果您为此编写自定义floodfill，则可以直接从中获取如下信息:

[pixels]中的卷

[pixels]中的表面积

边界框

位置(相对于肺)

重心或质心

所有这些都可以用作您的特征变量，也可以帮助进行拟合。

拟合找到的表面点

有很多方法可以解决这个问题，但是我会尽可能地简化它以提高性能和准确性。例如，您可以将质心用作椭球中心。然后从中找到min和max的远点，并将它们用作半轴(+/-一些正交校正)。然后只搜索这些初始值。有关更多信息，请参见:

How approximation search works

您可以在此处找到链接的 Q / A 中的使用示例。
[说明]

所有项目符号都适用于 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;