как извлечь границы изображения (изображение сканирования OCT / сетчатки)

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

оригинал

pic с границами

красная линия: верхняя граница первого слоя, зеленая линия: самая яркая линия во втором слое

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

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

    3 Solutions collect form web for “как извлечь границы изображения (изображение сканирования OCT / сетчатки)”

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

    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 или обрабатывать как подписанные значения … В моем примере я использовал предвзятость, так что серая область – это нулевой вывод … черный является самым негативным, а белый является наиболее положительным

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

    4. улучшить динамический диапазон

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

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

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

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

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

    предварительный просмотр

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

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

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

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

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

    • Как эффективно найти линию горизонта на высотной фотографии?
    • Обнаружение переломов в руке с использованием обработки изображений

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

    В коде используется Borland VCL с инкапсулированным GDI растровым изображением / холстом внизу (не важно, что вы просто визуализируете фактические настройки treshold) и свой собственный класс picture поэтому некоторые описания членов в порядке:

    • xs,ys разрешение
    • color p[ys][xs] прямой доступ к пикселям (32-битный формат пикселей, 8 бит на канал)
    • pf фактический выбранный формат пикселя для изображения см. ниже enum
    • enc_color/dec_color просто упаковывает распакованные цветовые каналы для разделения массива для удобной обработки нескольких пикселов … поэтому мне не нужно писать каждую функцию для каждого пиксельного формата отдельно
    • 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]=ci; } 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]=ci; } 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 { ci=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 { ci=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]; } //--------------------------------------------------------------------------- 

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

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

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

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

       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

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

    Полная рабочая реализация в 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 введите описание изображения здесь

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

    введите описание изображения здесь

     
    Interesting Posts for Van-Lav

    Что предпочтительнее использовать в Python: лямбда-функции или вложенные функции ('def')?

    заменяя только отдельные экземпляры символа регулярным выражением python

    Доступ к общедоступным страницам LinkedIn с использованием Python

    Shifted colorbar matplotlib

    Очень медленный поиск регулярных выражений

    Распараллеливать этот вложенный цикл for в python

    Найти имя переменной Python, которая была передана функции

    Отображение фотографий в ярлыках с помощью инструкции for с помощью tkinter, можно ли это сделать?

    Запросы Python, как ограничить полученный размер, скорость передачи и / или общее время?

    Элегантная питоническая кумса

    Как проверить, получены ли все данные с помощью Socket Socket в Python

    Как фильтровать перекрытие двух строк в большом файле в python

    Служба Windows Python – не отвечает на запуск / остановку из встроенного exe (но работает на python)

    Ошибка при загрузке DLL при импорте cv2

    Как я могу использовать результат возврата скрипта python PHP в реальном времени?

    Python - лучший язык программирования в мире.