matlab 上 3D 数据的椭球拟合

Question

我正在处理 3D 体积 的 CT 肺部图像，为了检测结节，我需要拟合椭圆体模型每个疑似结节，我怎样才能为那个做一个代码？？？结节是疑似肿瘤的对象，我的算法需要对每一个对象进行检查，将其近似为一个椭球体，从椭球体参数中计算出8个特征，构建一个分类器，通过训练和测试来判断是否为结节数据，所以我需要拟合这样的椭圆体

这里是体积的切片 CT 肺部图像

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

Answer 1

因为我们没有 3D 数据集，所以我从 2D.

开始

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

Fracture detection in hand using image proccessing

您可以在哪里找到您需要的所有解释。现在为了你的任务我会：

将RGBA转为灰度<0,765>

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

只需从边界的每条外线的中间向中间投射光线（扫描线），如果命中非白色像素则停止。我用 700 而不是 765 来补偿图像中的噪声。现在你得到了可用图像的边界框，所以裁剪掉其余部分。
计算直方图

为了补偿图像中的失真，使直方图足够平滑以消除所有不需要的噪声和间隙。然后从左向右求局部最大值（red）。这将用于二值化阈值（中间绿色）这是我的最终结果：
二值化图像

仅针对直方图中的**绿色*强度设置阈值图像。因此，如果 i0,i1 是直方图中左侧和右侧的局部最大强度，则阈值针对 (i0+i1)/2。这是结果：
除去肺以外的所有东西

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

应该没有留下黑色，剩下的白色可能是结节。
处理所有剩余的白色像素

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

现在您可以计算每个结节的特征。如果您为此编写自定义 floodfill 代码，那么您可以直接从中获取以下内容：
- 成交量 [pixels]
- 表面积 [pixels]
- 边界框
- 位置（相对于肺）
- 重心或质心
哪些都可以作为你的特征变量，也可以用来辅助拟合。
拟合找到的曲面点

有很多方法可以做到这一点，但我会尽可能简化它以提高性能和准确性。例如，您可以使用质心作为椭球中心。然后从它找到 min 和 max 远点并将它们用作半轴（+/- 一些正交校正）。然后搜索这些初始值。有关详细信息，请参阅：
您会在链接的 Q/A 中找到使用示例。

[备注]

所有子弹都适用于3D。在构造自定义 floodfill 时要小心递归尾部。过多的信息确实会很快 溢出您的堆栈 并且还会大大减慢速度。这里是我如何使用几个自定义 return 参数 + 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

您可能对我们为 open-source 软件 Icy http://icy.bioimageanalysis.org/

开发的最新插件感兴趣

插件名称为 FitEllipsoid，它允许通过首先单击正交视图上的几个点来快速将椭圆体拟合到图像内容。此处提供教程：https://www.youtube.com/watch?v=MjotgTZi6RQ

另请注意，我们在 GitHub 上提供了 Matlab 和 Java 源代码（但我无法提供它们，因为这是我第一次出现在该网站上）。

matlab 上 3D 数据的椭球拟合

Ellipsoid fitting for 3D data on matlab

3d

matlab

image-processing

dicom