基于matlab对图片中的函数图像重建

2017-06-25 本文已影响2333人未末星空

如需转载请先联系我

写在前面：
全部程序代码基于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