2016-06-21 5 views
3

У меня есть изображение (OCT), как показано ниже (оригинал). Как вы можете видеть, он имеет в основном 2 слоя. Я хочу создать изображение (показано на 3-м снимке), в котором красная линия указывает верхнюю границу 1-го слоя, зеленая показывает самую яркую часть 2-го слоя.как извлечь границы изображения (изображение сканирования OCT/сетчатки)

original

pic with borders

red line:top border of first layer, green line: brightest line in the 2nd layer

Я пытался просто пороги изображения. Затем я могу найти ребра, как показано на втором изображении. Но как можно производить красные/зеленые линии с этих границ?

PS: Я использую matlab (или OpenCV). Но любые идеи с любыми кодами языков/psudo приветствуются. спасибо заранее

+0

Какой язык (языки) вы используете? Добавьте соответствующие языковые теги, чтобы больше людей могли вам помочь. – rayryeng

+0

Я добавил информацию в исходное сообщение. Любые языки/идеи приветствуются. – wildcolor

ответ

4

Не слишком много времени на это прямо сейчас, но вы можете начать с этого:

  1. размыть изображение немного (удаление шума)

    простой свертка будет делать например, несколько раз с матрицей

    0.0 0.1 0.0 
    0.1 0.6 0.1 
    0.0 0.1 0.0 
    
  2. сделать вывод цвета по оси у

    производный цвет пикселя вдоль оси y ... например, я использовал это для каждого пикселя входного изображения:

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

    остерегайтесь результат подписан, так что вы можете нормализовать с помощью какой-то предвзятости или использовать abs значение или обрабатывать как знаковые ценности ... В моем примере я использовал уклон, так что серая область нулевого вывода ... черный является самым отрицательным и белым самым позитивным

  3. размытия изображения немного (удаление шума)

  4. повысить динамический диапазон

    Просто найти в изображении минимальное цвет c0 и максимальный цвет c1 и масштабировать все пиксели заданного диапазона <low,high>. Это сделает tresholding гораздо более стабильным на разных изображениях ...

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

    так, например, вы можете использовать low=0 и high=255

  5. Treshold все пиксели, которые больше, то Treshold

Полученное изображение выглядит так:

preview

Теперь просто:

  1. сегментировать красные участки
  2. удалить слишком маленькие участки
  3. термоусадочная/Перекрасить области оставить только среднюю точку на каждом x координат

    Верхняя часть точки красная, а нижняя - зеленая.

Это должно привести вас к очень близкому состоянию к вашему желаемому решению. Осторожно, размытие и деривация могут немного перемещать обнаруженные позиции с их реальной позиции.

Кроме того, для большего количества идей взглянуть на связанных с QA с:

[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 растровых/холст на день (не важен для вас только делают фактические параметры Treshold) и мой собственный picture класса, поэтому некоторые описания членов в порядке:

  • xs,ys разрешение
  • color p[ys][xs] прямого доступа пикселей (32-битный формат пикселей, так 8 бит на канал)
  • pf фактический формат выбран пиксель для изображения см ниже enum
  • enc_color/dec_color просто распакуйте цветные каналы в отдельный массив для удобной обработки нескольких пикселей ...поэтому мне не нужно писать каждую функцию для каждого PixelFormat отдельно
  • clear(DWORD c) заполняет изображение цвета c

color просто union из DWORD dd и BYTE db[4] и int i для простого доступа к каналу и или подписанных значений обработки.

Некоторые куски кода из него:

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

Я знаю его совсем немного кода ... и уравнения все, что вам нужно, но вы хотели это :) по своему усмотрению. Надеюсь, я не забыл что-то скопировать.

+0

Фантастический ответ! Также очень рад, что вы действительно продолжали прилагать усилия, чтобы ограничить свои очки! : D +1 для ya! –

+1

@AnderBiguri heh еще несколько парных ответов остается маркированным ... :) – Spektre

+0

@Spektre Не могли бы вы помочь объяснить (если возможно, какой-то код), как сделать вывод цвета по оси y? Делаю ли я производные вдоль оси y? Так что-то вроде этого (y2-y1)/(x2-x1)? Извините за глупый вопрос – wildcolor

3

Следующий набор инструкций (с использованием Matlab Image Processing Toolbox), кажется, работает хорошо для вашего изображения:

  1. Размытие изображение (Im) с медианным фильтром для уменьшения шума:

    ImB=medfilt2(Im,[20 20]); 
    
  2. Найти края с помощью Собела маски и растянуть их немного, чтобы соединить компоненты, близкие, и «чистые» общий имидж, чтобы избавиться от небольших участков

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

    EdgeClean.png

  3. затем У вас есть две зоны, которые можно обнаружить отдельно с помощью bwlabel

    L=bwlabel(EdgeClean); 
    
  4. Наконец, участок две зоны, соответствующие L = 1 и L = 2

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

    ImZones.png

Существует не так много шагов, потому что ваше начальное изображение довольно приятно, за исключением шума. Возможно, вам придется немного поиграть с параметрами, поскольку я начал с загруженной версии вашего изображения, которая может быть менее качественной, чем оригинальная. Конечно, это минимальный код, но я все еще надеюсь, что это поможет.

+0

Большое спасибо за фантастические ответы от «Спектр» и «Дpeурич». Оба они действительно помогают. Я новичок в обработке изображений. Благодарим за добавление кода Dpeurich. Я попробовал это на других моих изображениях. Кажется, что работает – wildcolor

+0

@wildcolor, если ответы на них вам помогут, вы должны принять его как решение своей проблемы (контрольная сумма около подсчета голосов). Также автором этого ответа является 'Toghe', а не' Dpeurich', или мне не хватает некоторых удаленных комментариев? – Spektre

+0

Spektre Спасибо за совет. Приносим извинения за неправильное имя @Toghe. Не знаю, почему я набрал это имя. – wildcolor

2

Полная рабочая реализация в октаве:

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

Идея заключается в том, чтобы обработать каждый столбец изображения в виде кривой (серых уровней) и ищет два пиков, каждый пик находится на границе слой.

Входное изображение оригинала связаны ОП: http://i.stack.imgur.com/PBnOj.pngenter image description here

изображения, сохраняемого в коде, как "dst.png" является следующее:

enter image description here