如何提取图像的边框（OCT/retinal 扫描图像）
how to extract the borders of an image (OCT/retinal scan image)
我有一个 (OCT) 图像，如下所示（原始）。如您所见，它主要有 2 层。我想生成一个图像（如图3所示），其中红线表示第一层的顶部边界，绿色表示第二层最亮的部分。
我试图简单地对图像进行阈值处理。然后我可以找到如图 2 所示的边缘。但是如何从这些边界产生 red/green 行呢？
PS：我正在使用 matlab（或 OpenCV）。但是欢迎使用任何 languages/psudo 代码的任何想法。提前致谢
现在没有太多时间，但你可以从这里开始：
稍微模糊图像（去除噪点）
简单的卷积可以用矩阵做几次
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
阈值大于阈值的所有像素
生成的图像是这样的：
现在只需：
- 分割红色区域
- 删除太小的区域
shrink/recolor 区域只留下每个 x 坐标的中点
最上面的点是红色，最下面的点是绿色。
这应该会让您非常接近您想要的解决方案。当心模糊和推导可以将检测到的位置从它们的真实位置移动一点。
还有更多想法请查看相关的 QAs:
- How to find horizon line efficiently in a high-altitude photo?
- Fracture detection in hand using image proccessing
[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 bitmap/canvas （对你来说不重要只是呈现实际的阈值设置) 和我自己的 picture class 所以一些成员的描述是有序的：
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 mask 找到边缘并稍微扩大它们以连接紧密的组件，'clean' 整体图像以去除小区域
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')
步骤不多，因为除了噪点之外，您的初始图像非常好。您可能需要对参数进行一些调整，因为我是从您的图像的下载版本开始的，它的质量可能低于原始图像。当然这是一个最小的代码，但我仍然希望这会有所帮助。
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"的图片如下：
我有一个 (OCT) 图像，如下所示（原始）。如您所见，它主要有 2 层。我想生成一个图像（如图3所示），其中红线表示第一层的顶部边界，绿色表示第二层最亮的部分。
我试图简单地对图像进行阈值处理。然后我可以找到如图 2 所示的边缘。但是如何从这些边界产生 red/green 行呢？
PS：我正在使用 matlab（或 OpenCV）。但是欢迎使用任何 languages/psudo 代码的任何想法。提前致谢
现在没有太多时间，但你可以从这里开始：
稍微模糊图像（去除噪点）
简单的卷积可以用矩阵做几次
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阈值大于阈值的所有像素
生成的图像是这样的：
现在只需：
- 分割红色区域
- 删除太小的区域
shrink/recolor 区域只留下每个
x坐标的中点最上面的点是红色，最下面的点是绿色。
这应该会让您非常接近您想要的解决方案。当心模糊和推导可以将检测到的位置从它们的真实位置移动一点。
还有更多想法请查看相关的 QAs:
- How to find horizon line efficiently in a high-altitude photo?
- Fracture detection in hand using image proccessing
[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 bitmap/canvas （对你来说不重要只是呈现实际的阈值设置) 和我自己的 picture class 所以一些成员的描述是有序的：
xs,ys分辨率color p[ys][xs]直接像素访问（32 位像素格式，因此每个通道 8 位）pf实际选择的图像像素格式参见下面的enumenc_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 mask 找到边缘并稍微扩大它们以连接紧密的组件，'clean' 整体图像以去除小区域
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')
步骤不多，因为除了噪点之外，您的初始图像非常好。您可能需要对参数进行一些调整，因为我是从您的图像的下载版本开始的，它的质量可能低于原始图像。当然这是一个最小的代码，但我仍然希望这会有所帮助。
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"的图片如下：