python - 如何提取图像的边界(OCT/视网膜扫描图像)

Question

我有一个 (OCT) 图像，如下所示(原始)。如您所见，它主要有 2 层。我想生成一个图像(如第 3 张图所示)，其中红线表示第 1 层的顶部边界，绿色表示第 2 层最亮的部分。

我试图简单地对图像进行阈值处理。然后我可以找到如图 2 所示的边缘。但是如何从这些边界产生红/绿线呢？

PS:我正在使用 matlab(或 OpenCV)。但是欢迎使用任何语言/伪代码的任何想法。提前致谢

Answer 1

现在没有太多时间，但您可以从这里开始:

1. 稍微模糊图像(去除噪点)

   简单的卷积可以用矩阵做几次

   0.0 0.1 0.0
   0.1 0.6 0.1
   0.0 0.1 0.0
   

2. 通过y轴进行颜色推导

   沿y导出像素颜色轴...例如，我将其用于输入图像的每个像素:

   pixel(x,y)=pixel(x,y)-pixel(x,y-1)
   

   注意结果是有符号的，因此您可以通过一些偏差进行归一化或使用 abs value 或 handle as signed values ...在我的示例中，我使用了偏差，因此灰色区域是零推导 ...黑色最负，白色最正

3. 稍微模糊图像(去除噪点)

4. 增强动态范围

   只需在图像中找到最小颜色 c0和最大颜色 c1并将所有像素重新缩放到预定义范围 <low,high> .这将使跨不同图像的阈值处理更加稳定...

   pixel(x,y)=low + ((pixel(x,y)-c0)*(high-low)/(c1-c0)
   

   例如，您可以使用 low=0和 high=255

5. 阈值大于阈值的所有像素

生成的图像是这样的:

现在只是:

1. 分割红色区域
2. 移除太小的区域

3. 缩小/重新着色区域以仅保留每个区域的中点 x坐标

   最上面的点是红色，最下面的点是绿色。

这应该会让您非常接近您想要的解决方案。当心模糊和推导会使检测到的位置稍微偏离其实际位置。

另外，如需更多想法，请查看相关的QA:

• How to find horizon line efficiently in a high-altitude photo?
• Fracture detection in hand using image proccessing

[Edit1] 我的 C++ 代码

picture pic0,pic1,pic2;     // pic0 is your input image, pic1,pic2 are output images
int x,y;
int tr0=Form1->sb_treshold0->Position;  // treshold from scrollbar (=100)
// [prepare image]
pic1=pic0;                  // copy input image pic0 to output pic1 (upper)
pic1.pixel_format(_pf_s);   // convert to grayscale intensity <0,765> (handle as signed numbers)
pic2=pic0;                  // copy input image pic0 to output pic2 (lower)

pic1.smooth(5);             // blur 5 times
pic1.derivey();             // derive colros by y
pic1.smooth(5);             // blur 5 times
pic1.enhance_range();       // maximize range

for (x=0;x<pic1.xs;x++)     // loop all pixels
 for (y=0;y<pic1.ys;y++)
  if (pic1.p[y][x].i>=tr0)  // if treshold in pic1 condition met
   pic2.p[y][x].dd=0x00FF0000; //0x00RRGGBB then recolor pixel in pic2

pic1.pixel_format(_pf_rgba); // convert the derivation signed grayscale to RGBA (as biased grayscale)

// just render actual set treshold
pic2.bmp->Canvas->Brush->Style=bsClear;
pic2.bmp->Canvas->Font->Color=clYellow;
pic2.bmp->Canvas->TextOutA(5,5,AnsiString().sprintf("Treshold: %i",tr0));
pic2.bmp->Canvas->Brush->Style=bsSolid;

代码在底部使用 Borlands VCL 封装的 GDI 位图/ Canvas (对您来说并不重要，只是呈现实际的阈值设置)和我自己的 picture类所以一些成员描述是有序的:

• xs,ys决议
• color p[ys][xs]直接像素访问(32 位像素格式，每个 channel 8 位)
• pf实际选择的图像像素格式请参阅 enum吼叫
• enc_color/dec_color只需将颜色 channel 打包解包到单独的数组，以便轻松处理多像素格式......所以我不需要为每个像素格式分别编写每个函数
• clear(DWORD c)用颜色填充图像 c

color只是union的 DWORD dd和 BYTE db[4]和 int i用于简单的 channel 访问和/或签名值处理。

其中的一些代码块:

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 _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
    };
//---------------------------------------------------------------------------
void dec_color(int *p,color &c,int _pf)
    {
    p[0]=0;
    p[1]=0;
    p[2]=0;
    p[3]=0;
         if (_pf==_pf_rgba) // 32 bit RGBA
         {
         p[0]=c.db[0];
         p[1]=c.db[1];
         p[2]=c.db[2];
         p[3]=c.db[3];
         }
    else if (_pf==_pf_s   ) // 32 bit signed int
         {
         p[0]=c.i;
         }
    else if (_pf==_pf_u   ) // 32 bit unsigned int
         {
         p[0]=c.dd;
         }
    else if (_pf==_pf_ss  ) // 2x16 bit signed int
         {
         p[0]=c.ii[0];
         p[1]=c.ii[1];
         }
    else if (_pf==_pf_uu  ) // 2x16 bit unsigned int
         {
         p[0]=c.dw[0];
         p[1]=c.dw[1];
         }
    }
//---------------------------------------------------------------------------
void dec_color(double *p,color &c,int _pf)
    {
    p[0]=0.0;
    p[1]=0.0;
    p[2]=0.0;
    p[3]=0.0;
         if (_pf==_pf_rgba) // 32 bit RGBA
         {
         p[0]=c.db[0];
         p[1]=c.db[1];
         p[2]=c.db[2];
         p[3]=c.db[3];
         }
    else if (_pf==_pf_s   ) // 32 bit signed int
         {
         p[0]=c.i;
         }
    else if (_pf==_pf_u   ) // 32 bit unsigned int
         {
         p[0]=c.dd;
         }
    else if (_pf==_pf_ss  ) // 2x16 bit signed int
         {
         p[0]=c.ii[0];
         p[1]=c.ii[1];
         }
    else if (_pf==_pf_uu  ) // 2x16 bit unsigned int
         {
         p[0]=c.dw[0];
         p[1]=c.dw[1];
         }
    }
//---------------------------------------------------------------------------
void enc_color(int *p,color &c,int _pf)
    {
    c.dd=0;
         if (_pf==_pf_rgba) // 32 bit RGBA
         {
         c.db[0]=p[0];
         c.db[1]=p[1];
         c.db[2]=p[2];
         c.db[3]=p[3];
         }
    else if (_pf==_pf_s   ) // 32 bit signed int
         {
         c.i=p[0];
         }
    else if (_pf==_pf_u   ) // 32 bit unsigned int
         {
         c.dd=p[0];
         }
    else if (_pf==_pf_ss  ) // 2x16 bit signed int
         {
         c.ii[0]=p[0];
         c.ii[1]=p[1];
         }
    else if (_pf==_pf_uu  ) // 2x16 bit unsigned int
         {
         c.dw[0]=p[0];
         c.dw[1]=p[1];
         }
    }
//---------------------------------------------------------------------------
void enc_color(double *p,color &c,int _pf)
    {
    c.dd=0;
         if (_pf==_pf_rgba) // 32 bit RGBA
         {
         c.db[0]=p[0];
         c.db[1]=p[1];
         c.db[2]=p[2];
         c.db[3]=p[3];
         }
    else if (_pf==_pf_s   ) // 32 bit signed int
         {
         c.i=p[0];
         }
    else if (_pf==_pf_u   ) // 32 bit unsigned int
         {
         c.dd=p[0];
         }
    else if (_pf==_pf_ss  ) // 2x16 bit signed int
         {
         c.ii[0]=p[0];
         c.ii[1]=p[1];
         }
    else if (_pf==_pf_uu  ) // 2x16 bit unsigned int
         {
         c.dw[0]=p[0];
         c.dw[1]=p[1];
         }
    }
//---------------------------------------------------------------------------
void picture::smooth(int n)
    {
    color   *q0,*q1;
    int     x,y,i,c0[4],c1[4],c2[4];
    bool _signed;
    if ((xs<2)||(ys<2)) return;
    for (;n>0;n--)
        {
        #define loop_beg for (y=0;y<ys-1;y++){ q0=p[y]; q1=p[y+1]; for (x=0;x<xs-1;x++) { dec_color(c0,q0[x],pf); dec_color(c1,q0[x+1],pf); dec_color(c2,q1[x],pf);
        #define loop_end enc_color(c0,q0[x  ],pf); }}
        if (pf==_pf_rgba) loop_beg for (i=0;i<4;i++) { c0[i]=(c0[i]+c0[i]+c1[i]+c2[i])>>2; clamp_u8(c0[i]);  } loop_end
        if (pf==_pf_s   ) loop_beg                   { c0[0]=(c0[0]+c0[0]+c1[0]+c2[0])/ 4; clamp_s32(c0[0]); } loop_end
        if (pf==_pf_u   ) loop_beg                   { c0[0]=(c0[0]+c0[0]+c1[0]+c2[0])>>2; clamp_u32(c0[0]); } loop_end
        if (pf==_pf_ss  ) loop_beg for (i=0;i<2;i++) { c0[i]=(c0[i]+c0[i]+c1[i]+c2[i])/ 4; clamp_s16(c0[i]); } loop_end
        if (pf==_pf_uu  ) loop_beg for (i=0;i<2;i++) { c0[i]=(c0[i]+c0[i]+c1[i]+c2[i])>>2; clamp_u16(c0[i]); } loop_end
        #undef loop_beg
        #define loop_beg for (y=ys-1;y>0;y--){ q0=p[y]; q1=p[y-1]; for (x=xs-1;x>0;x--) { dec_color(c0,q0[x],pf); dec_color(c1,q0[x-1],pf); dec_color(c2,q1[x],pf);
        if (pf==_pf_rgba) loop_beg for (i=0;i<4;i++) { c0[i]=(c0[i]+c0[i]+c1[i]+c2[i])>>2; clamp_u8(c0[i]);  } loop_end
        if (pf==_pf_s   ) loop_beg                   { c0[0]=(c0[0]+c0[0]+c1[0]+c2[0])/ 4; clamp_s32(c0[0]); } loop_end
        if (pf==_pf_u   ) loop_beg                   { c0[0]=(c0[0]+c0[0]+c1[0]+c2[0])>>2; clamp_u32(c0[0]); } loop_end
        if (pf==_pf_ss  ) loop_beg for (i=0;i<2;i++) { c0[i]=(c0[i]+c0[i]+c1[i]+c2[i])/ 4; clamp_s16(c0[i]); } loop_end
        if (pf==_pf_uu  ) loop_beg for (i=0;i<2;i++) { c0[i]=(c0[i]+c0[i]+c1[i]+c2[i])>>2; clamp_u16(c0[i]); } loop_end
        #undef loop_beg
        #undef loop_end
        }
    }
//---------------------------------------------------------------------------
void picture::enhance_range()
    {
    int i,x,y,a0[4],min[4],max,n,c0,c1,q,c;
    if (xs<1) return;
    if (ys<1) return;

    n=0;    // dimensions to interpolate
    if (pf==_pf_s   ) { n=1; c0=0; c1=127*3; }
    if (pf==_pf_u   ) { n=1; c0=0; c1=255*3; }
    if (pf==_pf_ss  ) { n=2; c0=0; c1=32767; }
    if (pf==_pf_uu  ) { n=2; c0=0; c1=65535; }
    if (pf==_pf_rgba) { n=4; c0=0; c1=  255; }

    // find min,max
    dec_color(a0,p[0][0],pf);
    for (i=0;i<n;i++) min[i]=a0[i]; max=0;
    for (y=0;y<ys;y++)
     for (x=0;x<xs;x++)
        {
        dec_color(a0,p[y][x],pf);
        for (q=0,i=0;i<n;i++)
            {
            c=a0[i]; if (c<0) c=-c;
            if (min[i]>c) min[i]=c;
            if (max<c) max=c;
            }
        }
    // change dynamic range to max
    for (y=0;y<ys;y++)
     for (x=0;x<xs;x++)
        {
        dec_color(a0,p[y][x],pf);
        for (i=0;i<n;i++) a0[i]=c0+(((a0[i]-min[i])*(c1-c0))/(max-min[i]+1));
//      for (i=0;i<n;i++) if (a0[i]<c0) a0[i]=c0; // clamp if needed
//      for (i=0;i<n;i++) if (a0[i]>c1) a0[i]=c1; // clamp if needed
        enc_color(a0,p[y][x],pf);
        }
    }
//---------------------------------------------------------------------------
void picture::derivey()
    {
    int i,x,y,a0[4],a1[4];
    if (ys<2) return;
    for (y=0;y<ys-1;y++)
     for (x=0;x<xs;x++)
        {
        dec_color(a0,p[y  ][x],pf);
        dec_color(a1,p[y+1][x],pf);
        for (i=0;i<4;i++) a0[i]=a1[i]-a0[i];
        enc_color(a0,p[y][x],pf);
        }
    for (x=0;x<xs;x++) p[ys-1][x]=p[ys-2][x];
    }
//---------------------------------------------------------------------------

我知道它有很多代码……方程式就是您所需要的，但您自己想要这个:)。希望我没有忘记复制一些东西。