Matlab上3D数据的椭球拟合

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

这里是体积的一部分CT 肺部图像

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

由于我们没有3D我开始使用的数据集2D.

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

• 使用图像处理检测手部骨折 https://stackoverflow.com/a/37046721/2521214

您可以在其中找到所需的所有解释。现在对于你的任务我会：

1. 将 RGBA 转为灰度<0,765>

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

2. 剪掉白色边框

   只需从边界每条外线的中间向中间投射光线（扫描线），并在遇到非白色像素时停止。我与700代替765来补偿图像中的噪声。现在您获得了可用图像的边界框，因此裁剪掉其余部分。

3. 计算直方图

   为了补偿图像中的失真，请平滑直方图以消除所有不需要的噪声和间隙。然后从左到右找到局部最大值（red）。这将用于二值化阈值（它们之间的中间值）green）这是我的最终结果：

4. 二值化图像

   只需根据直方图的**颗粒*强度对图像进行阈值处理即可。因此，如果i0,i1是直方图中左右两侧的局部最大强度，然后是阈值(i0+i1)/2。这是结果：

5. 除去肺以外的所有东西

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

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

6. 处理所有剩余的白色像素

   因此，只需循环遍历图像，并在第一个白色像素命中时用预定义的结节颜色或不同的对象索引填充它以供以后使用。我还区分了表面（aqua）和内部（magenta）。这是结果：

   现在您可以计算每个结节的特征。如果您编写自定义代码floodfill为此，您可以直接从中获取以下内容：

   • 音量输入[pixels]
   • 表面积在[pixels]
   • 边界框
   • 位置（相对于肺）
   • 重心或质心

   这些都可以用作特征变量，也可以帮助拟合。

7. 拟合找到的曲面点

   有很多方法可以实现这一点，但我会尽可能简化它以提高性能和准确性。例如，您可以使用质心作为椭球中心。然后找到min and max远离它的点并将它们用作半轴（+/-一些正交性校正）。然后就搜索这些初始值。欲了解更多信息，请参阅：

   • 近似搜索的工作原理 https://stackoverflow.com/a/36163847/2521214

   您将在链接中找到使用示例Q/As there.

[Notes]

所有项目符号都适用于3D。在构建自定义floodfill小心递归尾部。太多的信息会非常快溢出你的堆栈并且也会大大减慢速度。这是我如何使用一些自定义返回参数来处理它的小例子+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);
    }
//---------------------------------------------------------------------------

和这里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;