Clever Geek Handbook

Feature determination method for measuring the length of bent tubes in an image

I have hundreds of DNA nanotube images from fluorescence microscopy experiments and I would like to measure the tube length distribution in an automated way using image processing. Here is an example of a microscope image:

example DNA nanotube image

I have tried several methods for retrieving objects using python and skimage. I have tried using Canny's edge detection which successfully creates a diagram of each nanotube, however I am not clear on how to go from these contours to a specific measure of length. After applying Canny edge detection, I tried using Hough's probabilistic transform to fit straight lines to curves, which would make it easier to measure length. As you can see in these results:

Canny edge detection and hough transform results

Line substitution

is incompatible, and several lines are created in parallel for the same pipe structure.

Does anyone know of a direct method for measuring the length of these pipes?

+3


source to share


1 answer


I would start like this:

  • binarization image
  • find each set of tube pixels

  • fill fills this position with the tube color

    use any fill algorithm with 8 neighbors, and while filling, also highlight duplicate pixels in some counter cnt

    .

    if the size of the area is cnt

    too small, recolor it to the background, otherwise consider its size cnt/average_tube_width

    in the histogram.

Here's a simple C ++ example :

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;

Result for your input image:

result

The yellow lines are the length distribution (axis x

is tube length and y

is probability)

I am using my own image class for images, so some members:


xs,ys

- size in pixels
p[y][x].dd

- the pixel in the position (x,y)

as a 32-bit integer type
clear(color)

cleans the entire image using the color


resize(xs,ys)

resize the image to a new permit
bmp

VCL encapsulated GDI bitmap to Canvas

access
pf

contains the actual pixel image format:

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

and the pixels are encoded like this:



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


Ranges:

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

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

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

I also use my dynamic list template, so:


List<double> xxx;

matches double xxx[];


xxx.add(5);

add 5

to the end of the list
xxx[7]

an accessor array element (safe)
xxx.dat[7]

accessor array element (unsafe but fast direct access)
xxx.num

- the actual used array size
xxx.reset()

clears the array and sets xxx.num=0


xxx.allocate(100)

preset space for 100

elements

Now filling A * is done like this:

// these are picture:: members to speed up recursive fillings
int   _floodfill_rn;                                            // anti Qaru recursions
List<int> _floodfill_xy;                                        // anti Qaru 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;
    }

If you want to improve your score based on size, then if you have a multiple of avg, you have crossed tubes. So you can either try to calculate how many there are and account for the avg size for the histogram rather than use full size padding or * and find the endpoints. If you find more than two endpoints, you can try to distinguish between the tubes.

A * fill

So, first use A * fill to find the local max, and then A * fill that position again (so you start filling from the end point). Then find all the local max and based on the avg size and the actual tube size and number of endpoints, you can determine how many tubes are grouped together and how many are interconnected. Then you can try to make all possible combinations between the end points, and the closest to the average tube size for each tube is the most "correct" one. This should improve the accuracy much more.

If you do not know the tube thickness avg, you can use the A * fill for non-intersecting tubes directly to get the length. So in the second padding (from the end point), when the padding stops, the last padded ID is the tube length in pixels.

+2


source






More articles:


All Articles