logo标签列表

如何提取图像（OCT / 视网膜扫描图像）的边界

pythonmatlabopencvimage-processingimage-segmentation
3
我有一张类似于下面原始图片的光学相干断层扫描图像（OCT）。如您所见，它主要有两层。我想生成一张图片（如第三张图片所示），其中红线表示第一层的顶部边界，绿线显示第二层最亮的部分。
原始图片如下：original 带有边框的图片如下：pic with borders 红线：第一层的顶部边界，绿线：第二层最亮的部分的图片如下：red line:top border of first layer, green line: brightest line in the 2nd layer 我尝试简单地对图像进行了阈值处理。然后我可以找到像第二张图片中所示的边缘。但是，如何从这些边缘产生红/绿色线条呢？
PS: 我正在使用Matlab（或OpenCV）。但是，任何语言/伪代码的想法都受到欢迎。谢谢。
你使用的是哪种编程语言？请添加适当的语言标签，以便更多人帮助你。 - rayryeng
我已经将信息添加到原始帖子中。欢迎任何语言/想法。 - wildcolor
3个回答
4

目前我没有太多时间，但你可以从以下内容开始：

  1. blur the image a bit (remove noise)

    simple convolution will do for example few times with matrix

    0.0 0.1 0.0
0.1 0.6 0.1
0.0 0.1 0.0

  2. do color derivation by y axis

    derivate pixel color along y axis... for example I used this for each pixel of input image:

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

    beware the result is signed so you can normalize by some bias or use abs value or handle as signed values ... In my example I used bias so the gray area is zero derivation ... black is most negative and white is most positive

  3. blur the image a bit (remove noise)

  4. enhance dynamic range

    Simply find in image the min color c0 and max color c1 and rescale all pixels to predefined range <low,high>. This will make the tresholding much more stable across different images...

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

    so for example you can use low=0 and high=255

  5. treshold all pixels which are bigger then treshold

生成的图像如下所示：

preview

现在只需要：
  1. 分割红色区域
  2. 删除太小的区域

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

    最顶部的点是红色的，底部是绿色的。

这将使您接近所需的解决方案。请注意，模糊和导数可能会使检测到的位置与其真实位置有些偏差。
此外，更多想法，请查看相关的问答： [编辑1] 我的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;

这段代码使用Borland的VCL封装的GDI位图/画布（对您来说不重要，只是渲染实际阈值设置），以及我自己的picture类，因此需要一些成员变量的描述：

  • xs，ys分辨率
  • color p [ys] [xs]直接像素访问（32位像素格式，每通道8位）
  • pf图像的实际选定像素格式，请参见下面的enum
  • enc_color/dec_color将颜色通道打包解包到单独的数组中，以便轻松处理多个像素格式...因此，我不需要为每个像素格式分别编写每个函数
  • clear(DWORD c)用颜色c填充图像

color只是DWORD dd和BYTE db [4]以及int i的union，用于简单的通道访问和有符号值处理。

以下是其中的一些代码块：

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

我知道这是相当多的代码...和方程式，但你想要自己完成这个。希望我没有忘记复制什么东西。
太棒了！也非常高兴您确实继续努力限制您的项目符号！:+1: 给你！ - Ander Biguri
1
@AnderBiguri 嘿，还剩几百个答案要点出来... :) - Spektre
@wildcolor 是的，你可以沿着 y 轴求导 像素颜色... 例如，我对输入图像的每个像素使用了以下公式：pixel(x,y)=pixel(x,y)-pixel(x,y-1) 注意结果是有符号的，所以你可以通过一些偏差来进行归一化... 或者使用 ab 的值，或者将其处理为有符号的值... 在我的例子中，我使用了偏差，因此灰色区域的导数为零... 黑色是最负的，白色是最正的。 - Spektre
1
@wildcolor 是的，这3个步骤只是额外的后处理，以改善结果（如果需要），因为我不知道您得到的输入/输出的精度/质量。去除太小的区域很容易。只需对每个红色区域计算像素数（通过洪水填充重新着色或其他方式），并根据大小（像素数、边界框长宽比或其他）将其作为结果保留或通过重新着色为黑色等颜色来删除... - Spektre
1
还有一种基于模糊的方法...创建所有红点（最大强度）的掩码，其余部分为黑色（类似于ROI），现在对掩码进行几次模糊处理（越多，较大的区域将被去除），现在根据某个值阈值化掩码或仅保留最大强度像素。现在比较原始掩码和模糊掩码，并删除没有设置像素的区域。 - Spektre
显示剩余6条评论
3
以下是使用Matlab图像处理工具箱的一组指令，它似乎适用于您的图像：

  1. Blur your image (Im) with a median filter to decrease the noise:

    ImB=medfilt2(Im,[20 20]);

  2. Find the edges using sobel mask and dilate them a little bit to connect close components, and 'clean' the overall image to get rid of small areas

    edges = edge(ImB,'sobel');    
se = strel('disk',1);
EnhancedEdges = imdilate(edges, se);    
EdgeClean = bwareaopen(EnhancedEdges,1e3);

    EdgeClean.png

  3. You then have your two zones, that you can detect separately using bwlabel

    L=bwlabel(EdgeClean);

  4. Finally, plot the two zones corresponding to L=1 and L=2

    [x1 y1] = find(L==1);
[x2 y2] = find(L==2);
plot(y1,x1,'g',y2,x2,'r')

    ImZones.png

由于您的原始图像质量很好，除了噪点之外，所以步骤不多。您可能需要在参数上进行调整，因为我是从下载版本的图像开始的，这可能比原始图像质量低。当然，这只是一个最小化的代码，但我仍然希望能够帮到您。

非常感谢'Spektre'和'Dpeurich'提供的精彩答案。他们两个都真的很有帮助。我对图像处理还很陌生。感谢Dpeurich附上的代码。我已经在我的其他图片上尝试了一下，看起来效果不错 - wildcolor
@wildcolor 如果这个答案对你有用，你应该接受它作为解决问题的方案（在投票计数旁边的检查器）。此外，这个答案的作者是 Toghe 而不是 Dpeurich，或者我错过了一些已删除的评论吗？ - Spektre
Spektre，感谢您的建议。我为使用错误的名称@Toghe道歉。不确定为什么我会输入那个名字。 - wildcolor
实际上，在@wildcolor的第一条评论之后，我改了我的名字，所以你没有犯任何错误：），我道歉。 - Toghe
2
Octave的完整工作实现：
pkg load image
pkg load signal

median_filter_size = 10;
min_vertical_distance_between_layers = 35;
min_bright_level = 100/255;

img_rgb = imread("oct.png");% it's http://i.stack.imgur.com/PBnOj.png
img = im2double(img_rgb(:,:,1));
img=medfilt2(img,[median_filter_size median_filter_size]);

function idx = find_max(img,col_idx,min_vertical_distance_between_layers,min_bright_level)
  c1=img(:,col_idx);

  [pks idx]=findpeaks(c1,"MinPeakDistance",min_vertical_distance_between_layers,"MinPeakHeight",min_bright_level);

  if ( rows(idx) < 2 )
    idx=[1;1];
    return
  endif

  % sort decreasing peak value
  A=[pks idx];
  A=sortrows(A,-1);

  % keep the two biggest peaks
  pks=A(1:2,1);
  idx=A(1:2,2);

  % sort by row index
  A=[pks idx];
  A=sortrows(A,2);

  pks=A(1:2,1);
  idx=A(1:2,2);
endfunction

layers=[];
idxs=1:1:columns(img);
for col_idx=idxs
  layers=[layers find_max(img,col_idx,min_vertical_distance_between_layers,min_bright_level)];
endfor
hold on
imshow(img)
plot(idxs,layers(1,:),'r.')
plot(idxs,layers(2,:),'g.')

my_range=1:columns(idxs);

for i = my_range
  x = idxs(i);
  y1 = layers(1,i);
  y2 = layers(2,i);
  if  y1 > rows(img_rgb) || y2 > rows(img_rgb) || x > columns(img_rgb) || y1==1 || y2==1
    continue
  endif
  img_rgb(y1,x,:) = [255 0 0];
  img_rgb(y2,x,:) = [0 255 0];
endfor 

imwrite(img_rgb,"dst.png")

该想法是将图像的每一列处理为一条曲线（灰度级），并寻找两个峰值，每个峰值都在一个层的边界上。
输入图像是OP链接的原始图像：http://i.stack.imgur.com/PBnOj.png enter image description here 代码保存的图像为"dst.png"，如下所示：

enter image description here

网页内容由stack overflow 提供, 点击上面的
可以查看英文原文，
原文链接