用于测量图像中弯管长度的特征检测技术

Question

我有数百张荧光显微镜实验中的DNA纳米管的图像，我想用图像处理的自动方式测量管的长度分布。这是一个显微镜图像的例子：

​

我尝试了一些使用python和skimage的特征提取方法。我试过使用Canny边缘检测，它成功地创建了每个纳米管的轮廓，但我不清楚如何从这些轮廓到长度的最终测量。在应用Canny边缘检测之后，我尝试使用概率Hough变换将直线拟合到曲线上，这将使长度测量变得简单。正如您在这些结果中所看到的：

​

​

线管件不一致，并且为同一管结构并行创建了多条线。

有没有人知道一种简单的方法来测量这些管的长度？

Answer 1

我会这样开始：

1. binarize image
2. find tube
3. flood每组像素用管颜色填充位置

使用具有8个邻居的任何填充算法，并且在填充期间还计算一些计数器cnt中的重新着色的像素。

如果区域大小cnt太小，则将其重新着色到背景，否则将其大小cnt/average_tube_width计算到直方图中。

下面是一个简单的C++示例：

picture pic0,pic1;
    // pic0 - source img
    // pic1 - output img
//                    0xAARRGGBB
const DWORD col_backg=0x00202020;   // gray
const DWORD col_tube =0x00FFFFFF;   // white
const DWORD col_done =0x0000A0FF;   // aqua
const DWORD col_noise=0x00000080;   // blue
const DWORD col_error=0x00FF0000;   // red  (too smal _hist value)
const DWORD col_hist =0x00FFFF00;   // yellow
const DWORD col_test =0x01000000;   // A* filling start color (must be bigger then all other colors used)
int x,y,xx,yy,i;
DWORD c;
const int _hist=256;    // max area size for histogram
int hist[_hist];        // histogram  

// copy source image to output
pic1=pic0;
pic1.enhance_range();               // maximize dynamic range <0,255>^3
pic1.pixel_format(_pf_u);           // convert to grayscale <0,765>
pic1.threshold(100,766,col_backg,col_tube); // threshold intensity to binarize image
pic1.pf=_pf_rgba;                   // set as RGBA (without conversion)

// clear histogram
for (i=0;i<_hist;i++) hist[i]=0;
// find all tubes
for (y=0;y<pic1.ys;y++)
 for (x=0;x<pic1.xs;x++)
  if (pic1.p[y][x].dd==col_tube)
    {
    pic1.Astarfill(x,y,col_test);   // fill count area (8 neighbors)
    if (pic1._floodfill_n>5)        // if valid size
        {
        c=col_done;                 // set recolor color to done
        // update histogram
        if (pic1._floodfill_n<_hist) hist[pic1._floodfill_n]++;
         else c=col_error;
        }
    else c=col_noise;
    // recolor filled bbox with c
    for (yy=pic1._floodfill_y0;yy<=pic1._floodfill_y1;yy++)
     for (xx=pic1._floodfill_x0;xx<=pic1._floodfill_x1;xx++)
      if (pic1.p[yy][xx].dd>=col_test)
       pic1.p[yy][xx].dd=c;
    }
// render histogram
for (x=0;x<_hist;x++)
 for (i=0,y=pic1.ys-1;(y>=0)&&(i<hist[x]<<2);y--,i++)
   pic1.p[y][x].dd=col_hist;

输入图像的结果为：

​

​

黄线是长度分布(x轴是管长度，y是概率)

我对图片使用自己的图片类，所以一些成员是：

xs,ys是以像素为单位的图像大小

p[y][x].dd是(x,y)位置的像素，为32位整数类型

clear(color)使用color清除整个图像

resize(xs,ys)将图像大小调整为新分辨率

bmp是VCL封装的GDI位图，支持Canvas访问

pf保存图像的实际像素格式：

enum _pixel_format_enum
    {
    _pf_none=0, // undefined
    _pf_rgba,   // 32 bit RGBA
    _pf_s,      // 32 bit signed int
    _pf_u,      // 32 bit unsigned int
    _pf_ss,     // 2x16 bit signed int
    _pf_uu,     // 2x16 bit unsigned int
    _pixel_format_enum_end
    };

color和像素的编码方式如下：

union color
    {
    DWORD dd; WORD dw[2]; byte db[4];
    int i; short int ii[2];
    color(){}; color(color& a){ *this=a; }; ~color(){}; color* operator = (const color *a) { dd=a->dd; return this; }; /*color* operator = (const color &a) { ...copy... return this; };*/
    };

这些频段包括：

enum{
    _x=0,   // dw
    _y=1,

    _b=0,   // db
    _g=1,
    _r=2,
    _a=3,

    _v=0,   // db
    _s=1,
    _h=2,
    };

我也使用我的动态列表模板，所以：

List<double> xxx;与double xxx[];相同

xxx.add(5);将5添加到列表末尾

xxx[7]访问数组元素(safe)

xxx.dat[7]访问数组元素(不安全但快速直接访问)

xxx.num是数组的实际使用大小

xxx.reset()清除阵列并设置xxx.num=0

xxx.allocate(100)为100项目预分配空间

现在A*填充是这样实现的：

// these are picture:: members to speed up recursive fillings
int   _floodfill_rn;                                            // anti stack overflow recursions
List<int> _floodfill_xy;                                        // anti stack overflow pendng recursions
int   _floodfill_a0[4];                                         // recursion filled color and fill color
color _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;

// here the filling I used
void picture::Astarfill(int x,int y,DWORD id)
    {
    _floodfill_c0=p[y][x];
    _floodfill_c1.dd=id;
    _floodfill_n=0;
    _floodfill_x0=x;
    _floodfill_y0=y;
    _floodfill_x1=x;
    _floodfill_y1=y;
    _floodfill_rn=0;
    _floodfill_xy.num=0;

    if ((x<0)||(x>=xs)||(y<0)||(y>=ys)) return;

    int i;
    List<int> xy0,xy1,*p0,*p1,*pp;
    // first point
    p0=&xy0;
    p1=&xy1;
    p0->num=0;
    p0->add(x); p0->add(y); p[y][x].dd=id; _floodfill_n++;
    for (;p0->num;)
        {
        p1->num=0; id++;
        for (i=0;i<p0->num;)
            {
            x=p0->dat[i]; i++;
            y=p0->dat[i]; i++;
            x--; if ((x>=0)&&(y>=0)&&(x<xs)&&(y<ys)&&(p[y][x].dd==_floodfill_c0.dd)){ p1->add(x); p1->add(y); p[y][x].dd=id; _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; }
            y--; if ((x>=0)&&(y>=0)&&(x<xs)&&(y<ys)&&(p[y][x].dd==_floodfill_c0.dd)){ p1->add(x); p1->add(y); p[y][x].dd=id; _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; }
            x++; if ((x>=0)&&(y>=0)&&(x<xs)&&(y<ys)&&(p[y][x].dd==_floodfill_c0.dd)){ p1->add(x); p1->add(y); p[y][x].dd=id; _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; }
            x++; if ((x>=0)&&(y>=0)&&(x<xs)&&(y<ys)&&(p[y][x].dd==_floodfill_c0.dd)){ p1->add(x); p1->add(y); p[y][x].dd=id; _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; }
            y++; if ((x>=0)&&(y>=0)&&(x<xs)&&(y<ys)&&(p[y][x].dd==_floodfill_c0.dd)){ p1->add(x); p1->add(y); p[y][x].dd=id; _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; }
            y++; if ((x>=0)&&(y>=0)&&(x<xs)&&(y<ys)&&(p[y][x].dd==_floodfill_c0.dd)){ p1->add(x); p1->add(y); p[y][x].dd=id; _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; }
            x--; if ((x>=0)&&(y>=0)&&(x<xs)&&(y<ys)&&(p[y][x].dd==_floodfill_c0.dd)){ p1->add(x); p1->add(y); p[y][x].dd=id; _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; }
            x--; if ((x>=0)&&(y>=0)&&(x<xs)&&(y<ys)&&(p[y][x].dd==_floodfill_c0.dd)){ p1->add(x); p1->add(y); p[y][x].dd=id; _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; }
            }
        pp=p0; p0=p1; p1=pp;
        }
    _floodfill_rn=id-1;
    }

如果你想根据大小改进计数，那么如果你得到平均大小的倍数，你就会得到相交的管子。因此，您可以尝试计算它们的数量，并将平均大小计算到直方图中，而不是使用完整大小或使用A*填充并定位端点。如果找到两个以上的端点，可以尝试区分管状体。

​

​

因此，首先使用A* fill找到局部最大值，然后从该位置再次使用A* fill (所以您从端点开始填充)。然后找到所有局部maxs，并根据管的平均大小、实际大小和端点数确定有多少管组合在一起，以及有多少管是相互连接的。然后，您可以尝试在端点之间进行所有可能的组合，最接近每个管的平均管子大小的组合是最“正确”的组合。这应该会大大提高精确度。

如果您不知道管子的平均厚度，可以直接使用A*填充不相交的管子来获得长度。因此，在第二次填充(从端点开始)中，当填充停止时，最后填充的ID是以像素为单位的管的长度。