我有一张（OCT）图像，如下所示（原始）。如您所见，它主要有2层。我想制作一张图片（如第三张图片所示），其中红线表示第一层的上边框，绿色表示第二层最亮的部分。
我试图简单地对图像进行阈值处理。然后我可以找到第二张图片中显示的边缘。但是如何从这些边界产生红/绿线呢？
PS：我正在使用 matlab（或 OpenCV）。但是欢迎任何语言/伪代码的想法。提前致谢
【问题讨论】：
标签： python matlab opencv image-processing image-segmentation
现在没有太多时间，但你可以从这个开始：
稍微模糊一下图像（去除噪点）
简单的卷积可以用矩阵做几次
0.0 0.1 0.0
0.1 0.6 0.1
0.0 0.1 0.0
y轴颜色推导
沿y 轴导出像素颜色...例如，我将其用于输入图像的每个像素：
pixel(x,y)=pixel(x,y)-pixel(x,y-1)
注意结果是有符号的，因此您可以通过一些偏差进行归一化或使用abs 值或作为有符号值处理...在我的示例中，我使用了偏差，因此灰色区域为零导数...黑色是最负的，并且白色是最积极的
稍微模糊一下图像（去除噪点）
增强动态范围
只需在图像中找到最小颜色c0 和最大颜色c1 并将所有像素重新缩放到预定义范围<low,high>。这将使跨不同图像的阈值更加稳定...
pixel(x,y)=low + ((pixel(x,y)-c0)*(high-low)/(c1-c0)
例如，您可以使用low=0 和high=255
阈值大于阈值的所有像素
生成的图像是这样的：
现在只是：
缩小/重新着色区域以仅保留每个 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 位图/画布（对您而言并不重要，只是呈现实际的阈值设置）和我自己的 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];
}
//---------------------------------------------------------------------------
我知道它有相当多的代码......而方程式就是你所需要的，但你自己想要这个:)。希望我没有忘记复制一些东西。
【讨论】：
以下指令集（使用 Matlab 图像处理工具箱）似乎适用于您的图像：
使用中值滤波器模糊您的图像 (Im) 以减少噪点：
ImB=medfilt2(Im,[20 20]);
使用 sobel 掩码找到边缘并稍微扩大它们以连接紧密的组件，并“清理”整个图像以去除小区域
edges = edge(ImB,'sobel');
se = strel('disk',1);
EnhancedEdges = imdilate(edges, se);
EdgeClean = bwareaopen(EnhancedEdges,1e3);
然后你有你的两个区域，你可以使用 bwlabel 分别检测
L=bwlabel(EdgeClean);
最后，画出L=1和L=2对应的两个区域
[x1 y1] = find(L==1);
[x2 y2] = find(L==2);
plot(y1,x1,'g',y2,x2,'r')
没有太多步骤，因为您的初始图像除了噪点外非常好。您可能需要对参数进行一些调整，因为我从您下载的图像版本开始，该版本的质量可能低于原始图像。当然这是一个最小的代码，但我仍然希望这会有所帮助。
【讨论】：
Toghe而不是Dpeurich或者我错过了一些已删除的cmets？
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"的图片如下：
【讨论】：