利用matlab从图片中提取曲线坐标数据

小灰灰 2022-12-01 01:18 182阅读 0赞

### 目录 ###

*   *  0.引言
     *  1.思路详解与分析
     *  2.MATLAB程序

![在这里插入图片描述][watermark_type_ZmFuZ3poZW5naGVpdGk_shadow_10_text_aHR0cHM6Ly9ibG9nLmNzZG4ubmV0L3dlaXhpbl80Mzk2NDk5Mw_size_16_color_FFFFFF_t_70_pic_center]

## 0.引言 ##

在读文献的时，经常遇到这样的情况：文章里提出的方法好有趣啊，好想拿文中用的数据来试试看看能不能得到相近的结果，可是文中只有根据原始数据绘制的曲线图，没有数据。如下图所示。

![在这里插入图片描述][watermark_type_ZmFuZ3poZW5naGVpdGk_shadow_10_text_aHR0cHM6Ly9ibG9nLmNzZG4ubmV0L3dlaXhpbl80Mzk2NDk5Mw_size_16_color_FFFFFF_t_70_pic_center 1]

此时，如果能从文中把这幅图截取下来，输入到一个函数中去，最后能返回从图片中提取到的曲线的坐标数据，岂不美哉。这便是本文的工作。

## 1.思路详解与分析 ##

**1.1准备待提取数据的曲线图片**

将待提取数据的曲线的图片(如.jpg格式图片)，利用 `imread` 输入到 `matlab` 中。

**1.2曲线图片预处理与数据转换**

曲线图片预处理步骤的主要工作包含如下：

（1）图像二值化

将输入图像进行二值化处理，但分割得到的结果并不全为数据，其中可能还包括坐标轴等干扰点需要去除。

（2）获取从图片像素到曲线坐标的定标数据

首先，通过ginput()手动从图片中提取到两个像素点，这两个点分别为曲线坐标框的左上角和右下角。

![在这里插入图片描述][watermark_type_ZmFuZ3poZW5naGVpdGk_shadow_10_text_aHR0cHM6Ly9ibG9nLmNzZG4ubmV0L3dlaXhpbl80Mzk2NDk5Mw_size_16_color_FFFFFF_t_70_pic_center 2]

此时,便获得了曲线在图片上的像素范围

\[x\_*index\_min, x\_index\_max*\] & \[*y\_index\_min, y\_index\_max*\]

然后，手动输入实际曲线的数据坐标范围 \[x\_min, x\_max\] & \[ymin, y\_max\_\]. 如下所示。

![在这里插入图片描述][2020082513582621.png]

此时，一方面得到了像素坐标，一方面得到了实际坐标。接下来便利用这对数据，将图片中全部的像素坐标转换到实际坐标。

最终，得到了由图片提取到的数据散点图，如下:

![在这里插入图片描述][watermark_type_ZmFuZ3poZW5naGVpdGk_shadow_10_text_aHR0cHM6Ly9ibG9nLmNzZG4ubmV0L3dlaXhpbl80Mzk2NDk5Mw_size_16_color_FFFFFF_t_70_pic_center 3]

可以看到，此时得到的结果，虽然曲线与所需要的相近，但曲线外的部分，如坐标轴框、坐标轴刻度等也被转换成了数据，还需要进一步的处理。

**1.3数据的进一步处理并得到最终曲线**

这一步的主要工作是在数据中除去曲线之外的部分（包括坐标轴框、坐标轴刻度等）；以及解决一个x数据对应多个y数据的情况。

显然，坐标轴在整幅数据中，均处于边界位置，因此，很容易想到的一种方法是，设定阈值，将距离边界较近的数据直接删除。这里，设定了两个阈值，一个用来限定x方向上的数据，一个用来限定y方向上的数据。比如设定：rate\_x=0.08; rate\_y=0.05; 意思是阈值设定为曲线最前端8%和最后段8%的数据，曲线最顶端5%和最底端5%的数据。

进一步的，对于提取到的数据图，大多数情况一个x会对应若干y，因为数据是由图片转化而来，而图片的分辨率有限，一个实际数据点会用多个像素来表示。解决此问题的中心思想是将同一个x对应的若干个y取均值，但不能直接求均值，因为还有很多y是噪声(如坐标轴线、由图片噪声带入的干扰点等)需要先去除，在第一个问题中，通过限定y的范围，已经在一定程度上除去了部分干扰，在此基础上，我们求取一个x对应的这组y值的均值mean与标准差std，当某些y值位于\[mean-std , mean+std\]之外，则认为这些y值波动太大，将它们删去。

到这里，我们就将数据的处理部分基本完成了，我们将处理后的数据再次绘制成曲线，便可以得到如下

![在这里插入图片描述][watermark_type_ZmFuZ3poZW5naGVpdGk_shadow_10_text_aHR0cHM6Ly9ibG9nLmNzZG4ubmV0L3dlaXhpbl80Mzk2NDk5Mw_size_16_color_FFFFFF_t_70_pic_center 4]

对比处理之前的数据，由于限定了范围，因此曲线图片中带来的坐标框等内容转化而来的数据已经被删去。

将需要提取坐标的曲线图片，和我们提取并处理后的数据绘制的曲线，放在一起对比如下：  
![在这里插入图片描述][watermark_type_ZmFuZ3poZW5naGVpdGk_shadow_10_text_aHR0cHM6Ly9ibG9nLmNzZG4ubmV0L3dlaXhpbl80Mzk2NDk5Mw_size_16_color_FFFFFF_t_70_pic_center 5]

可以看到，与原曲线图片相比，提取到的数据曲线相似度能达到较高要求。但进一步观察会发现，右图曲线较左图而言，高频分量有一定的减少（即右图曲线更平滑），原因在于数据处理时，对同一个x对应的这组y值进行了均值处理，则在图像上近似反应为均值滤波，从而使得提取到的数据绘制成的曲线的高频分量被滤除。

最后，将提取到的最终数据，保存起来如下：

![在这里插入图片描述][watermark_type_ZmFuZ3poZW5naGVpdGk_shadow_10_text_aHR0cHM6Ly9ibG9nLmNzZG4ubmV0L3dlaXhpbl80Mzk2NDk5Mw_size_16_color_FFFFFF_t_70_pic_center 6]

**1.4进一步的讨论——曲线拟合**

通过对图片中曲线的数据提取，可以得到数值上的答案，这会带来进一步的思考，即能否得到这些数据的解析表达式。很容易想到，利用最小二乘法来拟合这些数据，这便涉及到了曲线的拟合。（插值与拟合可以这么理解：对于数据点集，若均落在曲线上，则该曲线为插值曲线，否则为拟合曲线）

对于一些简单的曲线图片(如下)，可以考虑用泰勒级数来近似，即多项式拟合。

![在这里插入图片描述][watermark_type_ZmFuZ3poZW5naGVpdGk_shadow_10_text_aHR0cHM6Ly9ibG9nLmNzZG4ubmV0L3dlaXhpbl80Mzk2NDk5Mw_size_16_color_FFFFFF_t_70_pic_center 7]

数据提取并拟合的结果展示如下

![在这里插入图片描述][watermark_type_ZmFuZ3poZW5naGVpdGk_shadow_10_text_aHR0cHM6Ly9ibG9nLmNzZG4ubmV0L3dlaXhpbl80Mzk2NDk5Mw_size_16_color_FFFFFF_t_70_pic_center 8]

同时还能得到拟合多项式的系数

![在这里插入图片描述][20200825140200890.png]

从而得到该曲线的多项式(这里采用四阶多项式)表达式为:

![在这里插入图片描述][20200825140231105.png_pic_center]

理论上，泰勒级数可以分解任何函数，但实际上，多项式拟合的次数太高，会出现龙格库塔现象，即摆尾现象。因此，多项式拟合的阶数不易过高，一般低于5阶。对于本文最开始的那幅曲线而言，仅5阶的泰勒级数就显得力不从心了，因此，对于这种存在波动剧烈的函数，可以考虑用傅里叶级数进行拟合，或者如果能提供先验知识，可以直接用先验表达式进行拟合。

在MATLAB中，提供了cftool工具箱，其提供了拟合与插值的GUI，使用非常方便，直接在命令窗输入cftool即可调用，cftool界面如下所示，其具体使用方法不在此赘述。

![在这里插入图片描述][watermark_type_ZmFuZ3poZW5naGVpdGk_shadow_10_text_aHR0cHM6Ly9ibG9nLmNzZG4ubmV0L3dlaXhpbl80Mzk2NDk5Mw_size_16_color_FFFFFF_t_70_pic_center 9]

## 2.MATLAB程序 ##

MATLAB源代码如下所示，和以往的风格一样，提供了详细的注释

% //提取图片中的曲线数据
    clear,clc,close all
    %% //图片与曲线间的定标
    im=imread('tu1.jpg'); %//读入图片(替换成需要提取曲线的图片)
    im=rgb2gray(im); %//灰度变化
    thresh = graythresh(im); %//二值化阈值
    im=im2bw(im,thresh); %//二值化
    set(0,'defaultfigurecolor','w')
    imshow(im) %//显示图片
    [y,x]=find(im==0); %//找出图形中的“黑点”的坐标。该坐标是一维数据。
    y=max(y)-y; %//将屏幕坐标转换为右手系笛卡尔坐标
    y=fliplr(y); %//fliplr()——左右翻转数组
    plot(x,y,'r.','Markersize', 2);
    disp('请在Figrure中先后点击实际坐标框的两个顶点(左上点和右下点)，即A、B两点. ');
    [Xx,Yy]=ginput(2); %//Xx,Yy——指实际坐标框的两个顶点
    min_x=input('最小的x值');  %//输入x轴最小值
    max_x=input('最大的x值');  %//输入x轴最大值
    min_y=input('最小的y值');  %//输入y轴最小值
    max_y=input('最大的y值');  %//输入y轴最大值
    x=(x-Xx(1))*(max_x-min_x)/(Xx(2)-Xx(1))+min_x;
    y=(y-Yy(1))*(min_y-max_y)/(Yy(2)-Yy(1))+max_y;
    plot(x,y,'r.','Markersize', 2);
    axis([min_x,max_x,min_y,max_y])  %//根据输入设置坐标范围
    title('由原图片得到的未处理散点图')
    %% //将散点转换为可用的曲线
    %//需处理的问题与解决思路
    %//(1)散点图中可能一个x对应好几个y <---> 保留mean()-std()到mean()+std()之间的y值 并取平均处理
    %//(2)曲线的最前端和最后段干扰较大 <---> 去掉曲线整体的前(如5%)和后5%
    %//(3)曲线的最顶端和最底段干扰较大 <---> 去掉曲线整体的上10%和下10%
    
    %//参数预设
    rate_x=0.08;  %//曲线的最前端和最后段删除比例
    rate_y=0.05;  %//曲线的最顶端和最底段删除比例
    
    [x_uni,index_x_uni]=unique(x);  %//找出有多少个不同的x坐标
    
    x_uni(1:floor(length(x_uni)*rate_x))=[];  %//除去前rate_x(如5%)的x坐标
    x_uni(floor(length(x_uni)*(1-rate_x)):end)=[];  %//除去后rate_x的x坐标
    index_x_uni(1:floor(length(index_x_uni)*rate_x))=[];  %//除去前rate_x的x坐标
    index_x_uni(floor(length(index_x_uni)*(1-rate_x)):end)=[];  %//除去后rate_x的x坐标
    
    [mxu,~]=size(x_uni);
    [mx,~]=size(x);
    for ii=1:mxu
        if ii==mxu
            ytemp=y(index_x_uni(ii):mx);
        else
            ytemp=y(index_x_uni(ii):index_x_uni(ii+1));
        end
        % //删除方差过大的异常点
        threshold1=mean(ytemp)-std(ytemp);
        threshold2=mean(ytemp)+std(ytemp);
        ytemp(find(ytemp<threshold1))=[];  %//删除同一个x对应的一段y中的异常点
        ytemp(find(ytemp>threshold2))=[];
        % //删除距顶端和底端较近的点
        thresholdy=(max_y-min_y)*rate_y;  %//y坐标向阈值
        ytemp(find(ytemp>max_y-thresholdy))=[];  %//删除y轴向距离顶端与底端距离小于rate_y的坐标
        ytemp(find(ytemp<min_y+thresholdy))=[];
        % //剩下的y求均值
        y_uni(ii)=mean(ytemp);
    end
    % //此时很多x_uni点处对应的y_uni为空,即NAN,要进一步删去这些空点
    x_uni(find(isnan(y_uni)))=[];
    y_uni(find(isnan(y_uni)))=[];
    % //画图
    figure,plot(x_uni,y_uni),title('经处理后得到的扫描曲线')
    axis([min_x,max_x,min_y,max_y])  %//根据输入设置坐标范围
    % //将最终提取到的x与y数据保存
    curve_val(1,:)=x_uni';
    curve_val(2,:)=y_uni;
    %% //对提取出的数据进行拟合(按实际情况进行修改)
    [p,s]=polyfit(curve_val(1,:),curve_val(2,:),4);  %//多项式拟合(为避免龙格库塔,多项式拟合阶数不宜太高)
    [y_fit,DELTA]=polyval(p,x_uni,s);  %//求拟合后多项式在x_uni对应的y_fit值
    figure,plot(x_uni,y_fit),title('拟合后的曲线')
    axis([min_x,max_x,min_y,max_y])  %//根据输入设置坐标范围

--------------------

转自[知乎专栏][Link 1]

[watermark_type_ZmFuZ3poZW5naGVpdGk_shadow_10_text_aHR0cHM6Ly9ibG9nLmNzZG4ubmV0L3dlaXhpbl80Mzk2NDk5Mw_size_16_color_FFFFFF_t_70_pic_center]: /images/20221123/eefb0b71c35940c8bf4fad6758bdd130.png
[watermark_type_ZmFuZ3poZW5naGVpdGk_shadow_10_text_aHR0cHM6Ly9ibG9nLmNzZG4ubmV0L3dlaXhpbl80Mzk2NDk5Mw_size_16_color_FFFFFF_t_70_pic_center 1]: /images/20221123/c48f3133775642668fb55f9aee7ba99f.png
[watermark_type_ZmFuZ3poZW5naGVpdGk_shadow_10_text_aHR0cHM6Ly9ibG9nLmNzZG4ubmV0L3dlaXhpbl80Mzk2NDk5Mw_size_16_color_FFFFFF_t_70_pic_center 2]: /images/20221123/8e346d5421d64a5180ea70dee4bd1bf6.png
[2020082513582621.png]: /images/20221123/01823a6a4aae494daea74fa211cc743e.png
[watermark_type_ZmFuZ3poZW5naGVpdGk_shadow_10_text_aHR0cHM6Ly9ibG9nLmNzZG4ubmV0L3dlaXhpbl80Mzk2NDk5Mw_size_16_color_FFFFFF_t_70_pic_center 3]: /images/20221123/b48d1ed0a9f74c209eda7d08de5df24d.png
[watermark_type_ZmFuZ3poZW5naGVpdGk_shadow_10_text_aHR0cHM6Ly9ibG9nLmNzZG4ubmV0L3dlaXhpbl80Mzk2NDk5Mw_size_16_color_FFFFFF_t_70_pic_center 4]: /images/20221123/5b619390d6a340589c03e8f0fd4e8f2d.png
[watermark_type_ZmFuZ3poZW5naGVpdGk_shadow_10_text_aHR0cHM6Ly9ibG9nLmNzZG4ubmV0L3dlaXhpbl80Mzk2NDk5Mw_size_16_color_FFFFFF_t_70_pic_center 5]: /images/20221123/9e95e7083d134a04afef9b1673af4441.png
[watermark_type_ZmFuZ3poZW5naGVpdGk_shadow_10_text_aHR0cHM6Ly9ibG9nLmNzZG4ubmV0L3dlaXhpbl80Mzk2NDk5Mw_size_16_color_FFFFFF_t_70_pic_center 6]: /images/20221123/c16a3fa4acd344c389c290f1722ec90c.png
[watermark_type_ZmFuZ3poZW5naGVpdGk_shadow_10_text_aHR0cHM6Ly9ibG9nLmNzZG4ubmV0L3dlaXhpbl80Mzk2NDk5Mw_size_16_color_FFFFFF_t_70_pic_center 7]: /images/20221123/7d62c9f6d33042a8adc5bfedbb493103.png
[watermark_type_ZmFuZ3poZW5naGVpdGk_shadow_10_text_aHR0cHM6Ly9ibG9nLmNzZG4ubmV0L3dlaXhpbl80Mzk2NDk5Mw_size_16_color_FFFFFF_t_70_pic_center 8]: /images/20221123/58c1126d3a2a4482b74e5d08c9ef0ae3.png
[20200825140200890.png]: /images/20221123/3daf84417ad741a1bba362f428405660.png
[20200825140231105.png_pic_center]: /images/20221123/430d2ce731934c3d8538c85e6401080b.png
[watermark_type_ZmFuZ3poZW5naGVpdGk_shadow_10_text_aHR0cHM6Ly9ibG9nLmNzZG4ubmV0L3dlaXhpbl80Mzk2NDk5Mw_size_16_color_FFFFFF_t_70_pic_center 9]: /images/20221123/b0ed2ca88d2444ee98a6bf20477a5be3.png
[Link 1]: https://zhuanlan.zhihu.com/p/52112012