如需转载请先联系我
写在前面：
全部程序代码基于Matlab 2015b 不保证版本过低的Matlab能正常使用
如果matlab版本太低没有内置以下函数请下载以下文件放到工作目录
已知问题：
matlab 7.0以下的版本不支持“~”这个符号以用于代替无关参数，所以在7.0版本以下的版本中无法实现本文用到的函数功能，请至少使用7.0以上的版本
"identifier" expected, "~" found. Error in
https://cn.mathworks.com/matlabcentral/newsreader/view_thread/324057
hough.m
https://drop.me/MYwpjo
houghpeaks.m
https://drop.me/BZwbWo
houghlines.m
https://drop.me/aGwy3a
题目预览
如果已知函数表达式，就可以利用计算机轻易地画出函数图形。相反，如果给你一张图片，里面有函数的图形，能不能获取其中的数据，从而进一步构建出函数的表达式呢？请以附件中的图片为例，解决以下几个问题：
- 问题1：如何获取图片中曲线的坐标数据？请以01.jpg为例，运用你的方法获取数据，并重新绘制图形。
- 问题2：如果图片中包含有坐标系，如何去除？请针对02.jpg、03.jpg应用你的方法进行处理，并运用1）的方法获取坐标数据。
- 问题3：当图片中有多条曲线时，如何获取每条曲线的坐标数据？请针对04.jpg、05.jpg应用你的方法，给出处理的结果。对05.jpg中获得的数据进行重新标定，使得标定后的数据的大小与原图中显示的结果尽可能的相符。
- 问题4：如何利用获得的坐标数据，重建函数表达式？请给出06.jpg、07.jpg中各曲线的函数表达式。
（注：在应用您的方法之前，可以适当地对图片做预处理，包括手工的处理，但应该进行说明，并遵循尽量减少手工处理的原则）
图片附件
01.jpg
02.jpg
03.jpg
04.jpg
05.jpg
06.jpg
07.jpg
解题思路
- 总体思想是去掉曲线以外的部分，对曲线重构，并进行坐标变换以尽可能接近原图中的函数
- 就本题而言，曲线以外的部分包括字母、数字、箭头，坐标轴（直线）
- 曲线重构通过二值化取点做
- 坐标变换通过屏幕坐标与笛卡尔坐标关系进行
- 获取函数表达式通过多项式拟合进行，考虑到简便性就不做幂函数、指数函数、对数函数的拟合，方法都是类似的
- 针对每幅图都有单独的m文件解题，但主要通过调用写好的函数。
程序源代码
模块函数（供调用）
findandpoly.m
%--------------------------------------------------------
%本程序用于对二值化后（黑曲线白背景）的图像进行曲线提取并拟合
%参数2个，BW（图像二值化后的矩阵），辅助参数他，分3种选择：0、普通、-1
%1是对曲线进行多项式拟合（3、5、7次）并和原提取到的曲线相对比（即画在同一张图里）
%2是进行指定n次多项式拟合，并与原曲线对比，同时在图中标注拟合曲线的表达式
%3是只进行坐标变换
%其他：无论被处理图像有没有坐标轴，执行此程序时都需要手动依次选择两个点，
%第一个是(0,0)，第二个是(0,1)，以用于重建坐标系并便于求曲线表达式
%如果需要其他特殊参数，请自行更改相关部分
%--------------------------------------------------------
function [p1,x,y]=findandpoly(BW,pic,t)
p1=[];
a=imread(pic);
%xx=size(a,1);%x1为图像横向长度，实际上这里用不上
yy=size(a,2);%y1为纵向长度
imshow(a);
[x,y] = ginput(1);%选(0,0)
x1=x;y1=yy-y;
[x,y] = ginput(1);%选(0,1)（对06.jpg是（2*pi,0））
x2=x;y2=yy-y;
%现在(x1,y1),(x2,y2)即为提取的点的笛卡尔坐标系的坐标
k=(y2-y1);%k是屏幕坐标系与图像坐标系的倍数
%k=(x2-x1)./(2.*pi);%针对06.jpg的处理改进
clear x;clear y;
[x,y]=find(BW==0);%找到黑色的点，并生成坐标（矩阵）
temp=x;%这里是进行图像坐标与屏幕坐标之间的变换
x=y;x=(x-x1)./k;%切换到和图片相符的坐标系，后续求出的拟合即是图片中的表达式了
y=(yy-temp);y=(y-y1)./k;
if (t~=-1)
clear x1;clear x2;clear y1;clear y2;
% mu = [mean(x); std(x)];%手动实现 centering and scaling：
%x1 = (x - mu(1))/mu(2);%即以横坐标的均值为中心、以横坐标的标准差来做归一化因子，将数据变得相对集中。
if (t~=0)
p1 = polyfit(x,y,t);%其中, x, y为已知数据点向量, 分别表示横,纵坐标, t为拟合多项式的次数
figure
plot(x,y,'k.',x,polyval(p1,x),'r')
legend('原始数据',[num2str(t) '次拟合数据'],2)%2表示图例在右上角
ftext=num2str(p1(t+1));
for i=1:t
ftext=[ftext '+' '(' num2str(p1(t-i+1)) '*x.^' num2str(i) ')'];
end
%ftext(end)=[];
ftext=strcat('y=',ftext)
p1=ftext;
title(ftext);
else
p1 = polyfit(x,y,3);%其中, x, y为已知数据点向量, 分别表示横,纵坐标, 3为拟合多项式的次数
p2 = polyfit(x,y,5);
p3 = polyfit(x,y,7);
figure
plot(x,y,'k.',x,polyval(p1,x),'r',x,polyval(p2,x),'g',x,polyval(p3,x),'b')
legend('原始数据','3次拟合数据','5次拟合数据','7次拟合数据',2)
end
axis equal%使坐标轴单位长度一直，否则图像会看起来扭曲
end
ginputpro.m
function [ii]=ginputpro(I)
figure
imshow(I);
%end
p=size(I,3);
if p==3
r=I(:,:,1);
g=I(:,:,2);
b=I(:,:,3);
I1=rgb2gray(I);
else
I1=I;
end
hold on
[x,y,c]=ginput(1);
m(1)=x;
n(1)=y;
plot(x,y,'r');
k=2;
while(c==1)
[x1,y1,c1]=ginput(1);
if c1==1
m(k)=x1;
n(k)=y1;
plot(x,y,'r');
line([m(k-1) m(k)],[n(k-1) n(k)]);
k=k+1;
c=c1;
else
break
end
end
line([m(k-1) m(1)],[n(k-1) n(1)]);
BW = roipoly(I1,m,n);
if p==3
r1=double(r).*double(BW);
g1=double(g).*double(BW);
b1=double(b).*double(BW);
i1=cat(3,r1,g1,b1);
figure
imshow(mat2gray(i1));
ii=i1;
else
i2=double(I1).*double(BW);
figure
imshow(mat2gray(i2));ii=i2;
end
ii=im2uint8(mat2gray(ii));
houghdiy.m
%-------------------------------------------------------------
%通过霍夫变换寻找直线，并将直线部分去掉（即将像素点RGB值赋为255成为白色）
%本函数输入参数两个，一个BW（图像二值化后的矩阵）、一个辅助参数t（用于确定是否对被检查到的直线范围进行扩大处理）
%-------------------------------------------------------------
function [BW]=houghdiy(BW,t)
I=~BW;
[H,theta,rho] = hough(BW);
P = houghpeaks(H,2,'threshold',ceil(0.5*max(H(:))));
lines = houghlines(BW,theta,rho,P,'FillGap',10,'MinLength',80);
for k = 1: length(lines)
xy = [lines(k).point1; lines(k).point2];
if xy(1,1)==xy(2,1)
I(xy(1,2):xy(2,2),(xy(1,1)-2):(xy(1,1)+2))=255;
else
syms kk bb x1 y1 x2 y2
[kk,bb]=solve('kk*x1+bb=y1','kk*x2+bb=y2','kk','bb');
kk=subs(kk,{x1,x2,y1,y2},{xy(1,1),xy(2,1),xy(1,2),xy(2,2)});
bb=subs(bb,{x1,x2,y1,y2},{xy(1,1),xy(2,1),xy(1,2),xy(2,2)});
kk=double(kk);
bb=double(bb);
if kk<=1 && kk>=-1
xx=linspace(xy(1,1),xy(2,1),abs(xy(1,1)-xy(2,1))+1);
yy=floor(kk*xx+bb);
else
yy=linspace(xy(1,2),xy(2,2),abs(xy(1,2)-xy(2,2))+1);
xx=floor((yy-bb)/kk);
end
for i=1:size(xx,2)
I((yy(i)-2):(yy(i)+2),xx(i))=255;
if (t==1)%有的时候检测到的直线部分不均匀，需要扩大处理范围
I((yy(i)-1):(yy(i)+2),xx(i))=255;
I((yy(i)-2):(yy(i)+3),xx(i))=255;
I((yy(i)-1):(yy(i)+3),xx(i))=255;
I((yy(i)-3):(yy(i)+2),xx(i))=255;
I((yy(i)-2):(yy(i)+1),xx(i))=255;
I((yy(i)-3):(yy(i)+1),xx(i))=255;
end
end
end
end
figure,imshow(I);%显示函数图像
BW=I;
quse.m
%-----------------------------------------------------
%im为输入的图像RGB数据，后面的三个数为要去除的颜色的RGB值
%-----------------------------------------------------
function im=quse(im,a,b,c)
jd = [a b c]; % 蓝色
njd = [1 1 1]*255; % 白色
tol = 1; % 阈值
for i = 1 : size(im, 1)
for j = 1 : size(im, 2);
temp = double(im(i, j, :)/255);
temp = temp(:).';
if norm(temp-jd) < tol
im(i, j, :) = njd;
im(i+1, j, :) = njd;
im(i, j+1, :) = njd;
im(i+1, j+1, :) = njd;
im(i-1, j, :) = njd;
im(i, j-1, :) = njd;
im(i-1, j-1, :) = njd;
end
end
end
imshow(im, []);
每幅图的主函数
Q01.m
%-------------------------------------------
%对于01：二值化->平滑处理->重绘图形->平滑处理结果与原图对比
%-------------------------------------------
close all
clear
pic='01.jpg';
a=imread(pic);
BW=im2bw(a,0.4);%二值化
[nop,x,y]=findandpoly(BW,pic,-1);
figure;plot(x,y);%默认打开句柄1，画出图像
print(2,'-djpeg','01_01.jpg');
title('图形重绘')
y(ceil(length(x)*rand(2,1))) = 3;%两种平滑处理，并画图对比
yy1 = smooth(x,y,0.1,'loess');
%yy2 = smooth(x,y,0.1,'rloess');第二种这个效果不好，舍去
[xx,ind] = sort(x);
figure;plot(xx,yy1(ind),'r-')%只重绘图形
title('图形平滑处理')
print(3,'-djpeg','01_02.jpg');
figure;plot(xx,y(ind),'b.',xx,yy1(ind),'y-')%包括与原二值图对比
title('图形对比')
print(4,'-djpeg','01_03.jpg');
Q02.m
%--------------------------------
%对于02：二值化->拟合对比->选取最优拟合并输出表达式
%--------------------------------
close all
clear
pic='02.jpg';
p=imread(pic);%读取图像
%I=im2bw(p,0.8);%直接二值化
BW=~im2bw(p,0.8);%后续处理图像（去掉直线）时是对I操作，即对原图二值化的图去直线
figure, imshow(BW);%显示二值化处理结果
title('二值化');
print(1,'-djpeg','02_01.jpg');
BW=houghdiy(BW,0);%霍夫变换，去直线
title('去直线')
print(2,'-djpeg','02_02.jpg');
findandpoly(BW,pic,0);%拟合并对比
title('多项式拟合对比');
print(3,'-djpeg','02_03.jpg');
[p1,~,~]=findandpoly(BW,pic,5);%拟合
title({p1 '最优拟合'});
print(4,'-djpeg','02_04.jpg');
%输出最后最优拟合系数p1
Q03.m
%-----------------------------------------
%对于03：预处理不需要的部分->去除指定颜色红色->去直线->二值化->拟合对比
%-----------------------------------------
close all
clear
pic='03.jpg';
im=imread(pic);
figure,imshow(im);
for t=1:3%选三次，依次是Y和旁边箭头；原点0；X和旁边箭头；当然也可以增加次数手动去除A、B
[x,y] = ginput(2); %在图像上手动选取两个顶点，从左上往右下选，矩形
for i = y(1):y(2)
for j = x(1):x(2)
im(i,j,1)=255;
im(i,j,2)=255;
im(i,j,3)=255;
end
end
end
%---------------------------------
%以上为选区域，做成白色（等同于去掉颜色）
%----------------------------------------
%第一步将手动选择区域处理为白色后，下面进行去除指定色彩（红色）
im=quse(im,1,0,0); % 红色
title('图像预处理');
print('-djpeg','03_01.jpg');
I=im2bw(im,0.8);
BW=~im2bw(im,0.8);
figure, imshow(BW);
title('二值化');
print('-djpeg','03_02.jpg');
BW=houghdiy(BW,0);%霍夫变换，去直线
title('去直线')
print('-djpeg','03_03.jpg');
findandpoly(BW,pic,0);%再次选点(0,0)(0,1);拟合并对比
title('多项式拟合对比');
print('-djpeg','03_04.jpg');
Q04.m
%————————————————————————————————————————
%对于04:预处理不需要的部分->去除蓝色->腐蚀图像（尽量减少直线部分）
%->去直线->再去除少数残留部分->拟合对比->选取最优拟合并输出表达式
%对于黑线步骤一致，只需改变需要去除的RGB值
%————————————————————————————————————————
close all
clear
pic='04.jpg';
im=imread(pic);
figure,imshow(im);
for t=1:4%选三次，依次是Y和旁边箭头；原点0；X和旁边箭头；当然也可以增加次数手动去除A、B
[x,y] = ginput(2); %在图像上手动选取两个顶点，从左上往右下选，矩形
for i = round(y(1)):round(y(2))
for j = round(x(1)):round(x(2))
im(i,j,1)=255;
im(i,j,2)=255;
im(i,j,3)=255;
end
end
end
im=quse(im,0,0,1);%去蓝色
title('图像预处理');
I=~im2bw(im,0.8);
se = strel('disk',1);%腐蚀图像
I2 = imerode(I,se);
figure,imshow(I2)
title('二值化');
print('-djpeg','04_01.jpg');
BW=houghdiy(I2,0);%霍夫变换，去直线
title('去直线')
print('-djpeg','04_02.jpg');
figure,imshow(BW);%显示处理后的图像
[x,y] = ginput(2); %需要手动再去除一点东西
for i = round(y(1)):round(y(2))
for j = round(x(1)):round(x(2))
BW(i,j)=255;
end
end
figure,imshow(BW);
title('进一步去直线')
print('-djpeg','04_03.jpg');
findandpoly(BW,pic,0);%选点(0,0)(0,1);拟合并对比
title('多项式拟合对比');
print('-djpeg','04_04.jpg');
[p1,~,~]=findandpoly(BW,pic,5);%指定拟合
title({p1 '黑线最优拟合'});
print('-djpeg','04_05.jpg');
%---------------以下为针对蓝线处理-----------------
im=imread(pic);
im=quse(im,0,0,0);%去黑色
title('图像预处理');
I=~im2bw(im,0.8);
se = strel('disk',1);%腐蚀图像
I2 = imerode(I,se);
figure,imshow(I2)
title('二值化');
print('-djpeg','04_06.jpg');
findandpoly(~I2,pic,0);%选点(0,0)(0,1);拟合并对比
title('多项式拟合对比');
print('-djpeg','04_07.jpg');
[p1,~,~]=findandpoly(~I2,pic,5);%指定拟合
title({p1 '蓝线最优拟合'});
print('-djpeg','04_08.jpg');
Q05.m
%-------------------------------------------------
%对于05：提取指定颜色+坐标变换，然后拟合即可
%这里以红线为例，提取红线更改输入的RGB值组合即可
%提取指定颜色->膨胀处理->二值化->拟合对比->选取最优拟合并输出表达式
%----------------------------------------------
close all
clear
pic='05.jpg';
a=imread(pic);
xx=size(a,1);%x1为图像横向长度
yy=size(a,2);%y1为纵向长度
imshow(a);
[x,y] = ginput(1);%选两个点即可，用于坐标变换
x1=x;y1=yy-y;
[x,y] = ginput(1);
x2=x;y2=yy-y;
k=(y2-y1);
a=double(a);
modeblue=zeros(size(a(:,:,1)));
modeblue=(a(:,:,1)>125).*(a(:,:,2)<125).*(a(:,:,3)<125);%<<>蓝色 ><<红色 <><黄色
I=~modeblue;
figure,imshow(I)
title('红线提取');
print('-djpeg','05_04.jpg');
se = strel('disk',1);%膨胀图像
I2 = imdilate(~I,se);
figure,imshow(I2);
title('二值化');
print('-djpeg','05_05.jpg');
[p1,~,~]=findandpoly(~I2,pic,9);%指定拟合
title({p1 '红线最优拟合'});
print('-djpeg','05_06.jpg');
Q06.m
%--------------------------------------
%对于06：
%预处理不需要的部分->去直线->拟合对比->选取最优拟合并输出表达式
%--------------------------------------
close all
clear
pic='06.jpg';
im=imread(pic);
figure,imshow(im);
for t=1:4%选四次，依次是Y和箭头；原点0；2pi;
[x,y] = ginput(2); %在图像上手动选取两个顶点，从左上往右下选，矩形
for i = round(y(1)):round(y(2))
for j = round(x(1)):round(x(2))
im(i,j,1)=255;
im(i,j,2)=255;
im(i,j,3)=255;
end
end
end
BW=~im2bw(im,0.8);
title('二值化');
print('-djpeg','06_01.jpg');
BW=houghdiy(BW,0);%霍夫变换，去直线
title('去直线')
print('-djpeg','06_02.jpg');
findandpoly(BW,pic,0);%选点(0,0)(2*pi,0);拟合并对比
title('多项式拟合对比');
print('-djpeg','06_03.jpg');
[p1,~,~]=findandpoly(BW,pic,9);%指定拟合
title({p1 '最优拟合'});
print('-djpeg','06_04.jpg');
Q07.m
%--------------------------------
%对于07：提取目标曲线->二值化->指定次数拟合
%--------------------------------
close all
clear
pic='07.jpg';
for kkk=1:7
kk=num2str(kkk);
im=imread(pic);
im=ginputpro(im);%多边形选取范围，并返回被选取的图像
title('选择要提取的目标曲线区域');
print('-djpeg',['07_' kk '_01.jpg']);
for i = 1 : size(im, 1)%把黑色部分填白，便于后续二值化及拟合
for j = 1 : size(im, 2);
if(im(i, j)==0)
im(i, j, :) = im(i, j, :)+255;
end
end
end
figure;imshow(im);
title('提取的曲线');
print('-djpeg',['07_' kk '_02.jpg']);
im=im2bw(im,0.8);%二值化
figure;imshow(im);
title('二值化并准备拟合');
print('-djpeg',['07_' kk '_03.jpg']);
[p1,~,~]=findandpoly(im,pic,3);%每一次都要重新选点，并指定拟合
title({p1 '指定为3次拟合'});
print('-djpeg',['07_' kk '_04.jpg']);
end
处理得到的图片
01
01_01.jpg
01_02.jpg
01_03.jpg
02
02_01.jpg
02_02.jpg
02_03.jpg
02_04.jpg
03
03_01.jpg
03_02.jpg
03_03.jpg
03_04.jpg
03_05.jpg
04
04_01.jpg
04_02.jpg
04_03.jpg
04_04.jpg
04_05.jpg
04_06.jpg
04_07.jpg
04_08.jpg
05
05_01.jpg
05_02.jpg
05_03.jpg
05_04.jpg
05_05.jpg
05_06.jpg
06
06_01.jpg
06_02.jpg
06_03.jpg
06_04.jpg
07
07_1_01.jpg
07_1_02.jpg
07_1_03.jpg
07_1_04.jpg
07_2_01.jpg
07_2_02.jpg
07_2_03.jpg
07_2_04.jpg
07_3_01.jpg
07_3_02.jpg
07_3_03.jpg
07_3_04.jpg
07_4_01.jpg
07_4_02.jpg
07_4_03.jpg
07_4_04.jpg
07_5_01.jpg
07_5_02.jpg
07_5_03.jpg
07_5_04.jpg
07_6_01.jpg
07_6_02.jpg
07_6_03.jpg
07_6_04.jpg
07_7_01.jpg
07_7_02.jpg
07_7_03.jpg
07_7_04.jpg
