how to extract the borders of an image (OCT/retinal scan image)

Question

I have an (OCT) image like shown below (original). As you can see, it mainly has 2 layers. I want to produce an image (shown in the 3rd picture), in which the red line indicates the top border of 1st layer, the green shows the brightest part of the 2nd layer.

I have tried to simply thresholded the image. Then I can find the edges like shown in the 2nd image. But how can produce the red/green lines from these borders?

PS: I am using matlab (or OpenCV). But any ideas with any languages/psudo codes are welcomed. thanks in advance

Answer 1

Do not have too much time for this right now but you can start with this:

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

The resulting image is like this:

Now just:

1. segment the red areas
2. remove too small areas

3. shrink/recolor areas to leave just the midpoint on each x coordinate

   The top most point is red and the bottom is green.

This should lead you to very close state to your wanted solution. Beware blurring and derivation can move the detected positions a bit from their real position.

Also for more ideas take a look at related QAs:

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

[Edit1] C++ code of mine for this

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;

The code is using Borlands VCL encapsulated GDI bitmap/canvas at bottom (not important for you just renders actual treshold settings) and my own picture class so some members description is in order:

• xs,ys resolution
• color p[ys][xs] direct pixel access (32bit pixel format so 8 bit per channel)
• pf actual selected pixel format for the image see the enum bellow
• enc_color/dec_color just pack unpack color channels to separate array for easy multi pixel format handling ... so I do not need to write each function for each pixelformat separately
• clear(DWORD c) fills image with color c

The color is just union of DWORD dd and BYTE db[4] and int i for simple channel access and or signed values handling.

Some chunks of code from it:

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

I know its quite a bit of code ... and the equations are all you need but you wanted this :) on your own. Hope I did not forget to copy something.

Answer 2

The following set of instructions (using Matlab Image processing toolbox) seems to work well for your image:

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);
   

   

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')
   

   

There are not a lot of steps because your initial image is quite nice except from the noise. You may need to play a bit on the parameters as I started from a downloaded version of your image which can be of lesser quality than the original one. Of course this is a minimal code, but I still hope this will help.

Answer 3

A full working implementation in 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 https://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")

The idea is to process every column of the image as a curve (of grey levels) and looking for two peaks, each peak is on the border of a layer.

The input image is the original linked by the OP: https://i.stack.imgur.com/PBnOj.png

The image saved by the code as "dst.png" is the following: