目的:由斑马鱼血管的二值图计算血管的交点。
方法:首先通过形态学计算血管骨架,然后手动选择血管骨架区域进行直线拟合,计算直线的斜率和截距,得到交点,该交点即为血管交点。
PS: 将直线叠加在原图上,结果可视化。
结果:
输入:
1血管二值图:
2 血管骨架:
3 血管二值+骨架:
血管骨架+拟合直线:
血管二值图+拟合直线:
视频原图+拟合直线:
代码:
img = imread('C:\Users\Admin\Desktop\1\gabor_wavelength_2.5重构+开运算.jpg');
% figure;imshow(img); title('openbw');
[h,w]=size(img);
for i =1:h
for j =1:w
if img(i,j)>200
bw(i,j)=0;
else
bw(i,j)=255;
end
end
end
for i =1:h
for j =1:w
if bw(i,j)==0
openbw(i,j)=0;
else
openbw(i,j)=1;
end
end
end
%% 去除面积小于1500的连通域
openbw=bwareaopen(openbw,1500,8);
% figure;imshow(openbw); title('openbw');
%% 填补黑色空洞
img_fill = imfill(openbw, 'holes');
% figure;imshow(img_fill); title('img_fill');
%% 开运算 平滑边缘
se = strel('disk',10);
openbw=imopen(img_fill,se);%对白色点而言
% figure;imshow(openbw); title('openbw');
%% MATLAB 自带接口函数 细化得到骨架
img_thining = bwmorph(openbw,'thin',Inf);
figure;imshow(img_thining);title('img_thining');
img_thining1=img_thining;
%% 细化的骨架叠加在二值图上
img_fill_thin=openbw;
for i =1:h
for j =1:w
if img_thining(i,j)==1
img_fill_thin(i,j)=0;
end
end
end
% figure;imshow(img_fill_thin); title('img_fill_thin');
%% 手动框选骨架区域,进行直线拟合,血管1
handle = imrect;
position = getPosition(handle);% 拖动鼠标获得兴趣区域,pos有四个值,兴趣区域左上角的像素坐标和区域的长宽
x1 = floor(position(1));
y1 = floor(position(2));
x2 = floor(position(3)+position(1));
y2 = floor(position(4)+position(2));
n=1;
%扫描每一个像素,并记录白点(值为1)坐标及个数
for i =y1:y2
for j =x1:x2
if img_thining(i,j)==1
points_for_line(n,1)=j;
points_for_line(n,2)=i;
n=n+1;
end
end
end
% figure;imshow(img_thining); title('img_thining——test');
t=polyfit(points_for_line(:,1),points_for_line(:,2),1);
%最小二乘法拟合直线
A = 0.0;
B = 0.0;
C = 0.0;
D = 0.0;
for i=1:1:n-1
A=A+points_for_line(i,1)*points_for_line(i,1);
B=B+points_for_line(i,1);
C=C+points_for_line(i,1)*points_for_line(i,2);
D=D+points_for_line(i,2);
end
%直线斜率和截距
a1 = (C*n - B*D) / (A*n - B*B);
b1 = (A*D - C*B) / (A*n - B*B);
%% 手动框选骨架区域,进行直线拟合,血管2
imshow(img_thining);
handle = imrect;
position = getPosition(handle);% 拖动鼠标获得兴趣区域,pos有四个值,兴趣区域左上角的像素坐标和区域的长宽
x1 = floor(position(1));
y1 = floor(position(2));
x2 = floor(position(3)+position(1));
y2 = floor(position(4)+position(2));
n=1;
%扫描每一个像素,并记录白点(值为1)坐标及个数
for i =y1:y2
for j =x1:x2
if img_thining(i,j)==1
points_for_line(n,1)=j;
points_for_line(n,2)=i;
n=n+1;
end
end
end
t=polyfit(points_for_line(:,1),points_for_line(:,2),1);
plot(points_for_line(:,1),points_for_line(:,2),'*',points_for_line(:,1),polyval(t,points_for_line(:,1)))
%最小二乘法拟合直线
A = 0.0;
B = 0.0;
C = 0.0;
D = 0.0;
for i=1:1:n-1
A=A+points_for_line(i,1)*points_for_line(i,1);
B=B+points_for_line(i,1);
C=C+points_for_line(i,1)*points_for_line(i,2);
D=D+points_for_line(i,2);
end
%直线斜率和截距
a2 = (C*n - B*D) / (A*n - B*B);
b2 = (A*D - C*B) / (A*n - B*B);
%% 将拟合直线叠加在骨架上
%二值图转rgb彩色图片
for i =1:h
for j =1:w
if img_thining1(i,j)==1
img_thining_rgb(i,j)=255;
else
img_thining_rgb(i,j)=0;
end
end
end
imrgb_thining=repmat(img_thining_rgb,[1,1,3]);
%画线1,线宽3,红色
for i=1:1:w
y=a1*i+b1;
y2=int32(y); %把y转换成整数
%把线上的通道1(红色)置为255
imrgb_thining(y2-1,i,1)=255;
imrgb_thining(y2,i,1)=255;
imrgb_thining(y2+1,i,1)=255;
%把线上的其他通道(绿色和蓝色)置为0
imrgb_thining(y2-1,i,2)=0;
imrgb_thining(y2,i,2)=0;
imrgb_thining(y2+1,i,2)=0;
imrgb_thining(y2-1,i,3)=0;
imrgb_thining(y2,i,3)=0;
imrgb_thining(y2+1,i,3)=0;
end
%画线2,线宽3,红色
for i=1:1:w
y=a2*i+b2;
y2=int32(y); %把y转换成整数
if y2>2 && y2<480
%把线上的通道1(红色)置为255
imrgb_thining(y2-1,i,1)=255;
imrgb_thining(y2,i,1)=255;
imrgb_thining(y2+1,i,1)=255;
%把线上的其他通道(绿色和蓝色)置为0
imrgb_thining(y2-1,i,2)=0;
imrgb_thining(y2,i,2)=0;
imrgb_thining(y2+1,i,2)=0;
imrgb_thining(y2-1,i,3)=0;
imrgb_thining(y2,i,3)=0;
imrgb_thining(y2+1,i,3)=0;
end
end
%显示图片
imshow(imrgb_thining);
%% 将拟合直线叠加在血管二值图上
%二值图转rgb彩色图片
for i =1:h
for j =1:w
if openbw(i,j)==1
img_fill_rgb(i,j)=255;
else
img_fill_rgb(i,j)=0;
end
end
end
imrgb_fill=repmat(img_fill_rgb,[1,1,3]);
%画线1,线宽3,红色
for i=1:1:w
y=a1*i+b1;
y2=int32(y); %把y转换成整数
%把线上的通道1(红色)置为255
imrgb_fill(y2-1,i,1)=255;
imrgb_fill(y2,i,1)=255;
imrgb_fill(y2+1,i,1)=255;
%把线上的其他通道(绿色和蓝色)置为0
imrgb_fill(y2-1,i,2)=0;
imrgb_fill(y2,i,2)=0;
imrgb_fill(y2+1,i,2)=0;
imrgb_fill(y2-1,i,3)=0;
imrgb_fill(y2,i,3)=0;
imrgb_fill(y2+1,i,3)=0;
end
%画线2,线宽3,红色
for i=1:1:w
y=a2*i+b2;
y2=int32(y); %把y转换成整数
if y2>2 && y2<480
%把线上的通道1(红色)置为255
imrgb_fill(y2-1,i,1)=255;
imrgb_fill(y2,i,1)=255;
imrgb_fill(y2+1,i,1)=255;
%把线上的其他通道(绿色和蓝色)置为0
imrgb_fill(y2-1,i,2)=0;
imrgb_fill(y2,i,2)=0;
imrgb_fill(y2+1,i,2)=0;
imrgb_fill(y2-1,i,3)=0;
imrgb_fill(y2,i,3)=0;
imrgb_fill(y2+1,i,3)=0;
end
end
%显示图片
imshow(imrgb_fill);
%% 将拟合直线叠加在视频原图上
img_origin=imread('C:\Users\Admin\Desktop\斑马鱼10.30\bandicam 2020-10-30 18-47-23-168_origin\001.jpg');
imrgb_origin=img_origin;
%画线1,线宽3,红色
for i=1:1:w
y=a1*i+b1;
y2=int32(y); %把y转换成整数
%把线上的通道1(红色)置为255
imrgb_origin(y2-1,i,1)=255;
imrgb_origin(y2,i,1)=255;
imrgb_origin(y2+1,i,1)=255;
%把线上的其他通道(绿色和蓝色)置为0
imrgb_origin(y2-1,i,2)=0;
imrgb_origin(y2,i,2)=0;
imrgb_origin(y2+1,i,2)=0;
imrgb_origin(y2-1,i,3)=0;
imrgb_origin(y2,i,3)=0;
imrgb_origin(y2+1,i,3)=0;
end
%画线2,线宽3,红色
for i=1:1:w
y=a2*i+b2;
y2=int32(y); %把y转换成整数
if y2>2 && y2<480
%把线上的通道1(红色)置为255
imrgb_origin(y2-1,i,1)=255;
imrgb_origin(y2,i,1)=255;
imrgb_origin(y2+1,i,1)=255;
%把线上的其他通道(绿色和蓝色)置为0
imrgb_origin(y2-1,i,2)=0;
imrgb_origin(y2,i,2)=0;
imrgb_origin(y2+1,i,2)=0;
imrgb_origin(y2-1,i,3)=0;
imrgb_origin(y2,i,3)=0;
imrgb_origin(y2+1,i,3)=0;
end
end
%显示图片
imshow(imrgb_origin);
%% 计算直线交点
x=(b2-b1)/(a1-a2);
y=a1*x+b1;
disp('x=');
disp(x);
disp('y=');
disp(y);
%% 保存到文件夹中
tmppath=pwd; %保存当前工作目录
cd('C:\Users\Admin\Desktop\1\'); %把当前工作目录切换到图片存储文件夹
fname='1血管二值图.jpg';
imwrite(openbw,fname);
fname='1血管骨架.jpg';
imwrite(img_thining,fname);
fname='1血管二值+骨架.jpg';
imwrite(img_fill_thin,fname);
fname='1血管骨架+拟合直线.jpg';
imwrite(imrgb_thining,fname);
fname='1血管二值图+拟合直线.jpg';
imwrite(imrgb_fill,fname);
fname='1视频原图+拟合直线.jpg';
imwrite(imrgb_origin,fname);
cd(tmppath) %切回原工作目录