找回密码
 注册
关于网站域名变更的通知
查看: 575|回复: 1
打印 上一主题 下一主题

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

[复制链接]

该用户从未签到

跳转到指定楼层
1#
发表于 2021-1-20 18:57 | 只看该作者 |只看大图 回帖奖励 |倒序浏览 |阅读模式

EDA365欢迎您登录!

您需要 登录 才可以下载或查看,没有帐号?注册

x

: M2 |" y% W1 I% ^- V目录
& [9 K3 i9 x+ ~( a0.引言
) q( ^* D4 C/ U7 ~, C+ n" I  l7 x. w1.思路详解与分析3 v9 {$ Y6 u/ K* U: N
2.MATLAB程序
' b$ p1 S& ]) X % X* Y& V/ X% {' M# @9 _9 @2 m  L
0.引言6 @7 K0 u5 h- v2 f
  在读文献的时,经常遇到这样的情况:文章里提出的方法好有趣啊,好想拿文中用的数据来试试看看能不能得到相近的结果,可是文中只有根据原始数据绘制的曲线图,没有数据。如下图所示。
/ O4 j/ }, l; V4 Q( J5 P2 D+ |' ^
& G0 Y$ \2 F& {/ b4 w4 ~ " Z7 X7 _+ D# k6 O& X

5 z, A4 I* h$ W  X/ M  j+ K  此时,如果能从文中把这幅图截取下来,输入到一个函数中去,最后能返回从图片中提取到的曲线的坐标数据,岂不美哉。这便是本文的工作。9 {4 c4 |( l) E2 `9 G: G7 ]
2 f& Y% i2 ?& H  B: j
1.思路详解与分析/ X$ N& M' K) a$ {  {" `- {1 P. y
1.1准备待提取数据的曲线图片
; L) g# O% u8 d' i- G& O( e3 l" N
5 `! }; I5 v) f4 c, W  x; _将待提取数据的曲线的图片(如.jpg格式图片),利用 imread 输入到 matlab 中。
) x) L( O+ ~; M. n) v/ s2 U, r$ U2 m$ \) I9 B
1.2曲线图片预处理与数据转换
6 y$ w0 _; F+ W' S5 `* i  a$ ]  v; J
曲线图片预处理步骤的主要工作包含如下:
& |+ d( z% }! k0 s. c. J/ V
' |: W  C$ E; ?0 \0 R(1)图像二值化
, Q4 }- a, f% l4 G- Q/ c
/ x. F3 @2 G& a, U将输入图像进行二值化处理,但分割得到的结果并不全为数据,其中可能还包括坐标轴等干扰点需要去除。: i3 e: ]5 e+ [) o, a1 I# m
- H# r1 r& G# K' `
(2)获取从图片像素到曲线坐标的定标数据8 U3 V3 D! ~9 o* a

" a+ b2 i, I+ h5 K9 R$ f6 ^0 M5 e! V首先,通过ginput()手动从图片中提取到两个像素点,这两个点分别为曲线坐标框的左上角和右下角。
& _4 g1 V% j& H# o5 W# \9 |8 V" L$ H
+ Q4 [0 Q# ?) J/ Y( k 8 i2 F; a  P2 y; M/ X; r

2 C1 v+ ]/ _$ @, T! L' x6 Q- p此时,便获得了曲线在图片上的像素范围+ b( a+ [2 K$ T1 _; Y4 V& X  z9 J
* W9 o; D' ?  b6 W! P/ e2 H  @
[x_index_min, x_index_max] & [y_index_min, y_index_max]
4 ]: C5 w. U# [0 Q# ]% v6 F, F7 X  C) Q
然后,手动输入实际曲线的数据坐标范围 [x_min, x_max] & [ymin, y_max_]. 如下所示。
0 S% F. A, O3 X8 B* T+ N( Q
- H  u& u; `9 Z' m1 D  ?5 u$ L
# c0 j+ Z, b+ d/ E4 F1 H2 w& @2 G$ ^& D& N' g5 q
此时,一方面得到了像素坐标,一方面得到了实际坐标。接下来便利用这对数据,将图片中全部的像素坐标转换到实际坐标。' T, K! v  V9 i3 V3 Y" y

0 F- G% L) J+ Y最终,得到了由图片提取到的数据散点图,如下:# P0 u0 V# \! z: H1 Y; o
" h1 o9 @6 [& X* E
$ l* P0 k1 I6 c* E7 l* C# T  K5 ^2 \
可以看到,此时得到的结果,虽然曲线与所需要的相近,但曲线外的部分,如坐标轴框、坐标轴刻度等也被转换成了数据,还需要进一步的处理。
8 Z) w7 {/ E( ]
( c/ u# A" U) Q' J1.3数据的进一步处理并得到最终曲线9 T) T( k3 E" r+ L  U3 B- A/ J, `

9 W3 G; g) o# Y# Q3 g- B- I这一步的主要工作是在数据中除去曲线之外的部分(包括坐标轴框、坐标轴刻度等);以及解决一个x数据对应多个y数据的情况。; K/ ]* O% I8 F" U/ ?: q
) [  `8 W1 ]( v4 }$ D
显然,坐标轴在整幅数据中,均处于边界位置,因此,很容易想到的一种方法是,设定阈值,将距离边界较近的数据直接删除。这里,设定了两个阈值,一个用来限定x方向上的数据,一个用来限定y方向上的数据。比如设定:rate_x=0.08; rate_y=0.05; 意思是阈值设定为曲线最前端8%和最后段8%的数据,曲线最顶端5%和最底端5%的数据。2 k4 t5 [4 p) r/ X, a
" w" E$ s7 Z2 J5 @+ w/ ^+ o
进一步的,对于提取到的数据图,大多数情况一个x会对应若干y,因为数据是由图片转化而来,而图片的分辨率有限,一个实际数据点会用多个像素来表示。解决此问题的中心思想是将同一个x对应的若干个y取均值,但不能直接求均值,因为还有很多y是噪声(如坐标轴线、由图片噪声带入的干扰点等)需要先去除,在第一个问题中,通过限定y的范围,已经在一定程度上除去了部分干扰,在此基础上,我们求取一个x对应的这组y值的均值mean与标准差std,当某些y值位于[mean-std , mean+std]之外,则认为这些y值波动太大,将它们删去。
6 q& e6 N4 l" _" e( h/ K
  D; X$ g" T2 D- C1 l3 B到这里,我们就将数据的处理部分基本完成了,我们将处理后的数据再次绘制成曲线,便可以得到如下
( Q0 z  x6 Z  T# I1 x( \$ `1 {4 u/ P' @1 {! M3 l5 f% ?5 L

! e& Y% [; B; I, Y: a9 R
8 Q- S6 s/ k& B' }0 c1 `对比处理之前的数据,由于限定了范围,因此曲线图片中带来的坐标框等内容转化而来的数据已经被删去。: x% o/ ]9 `2 R
# }$ w2 a$ m! r* r) l
将需要提取坐标的曲线图片,和我们提取并处理后的数据绘制的曲线,放在一起对比如下:' E" V3 t! N  R

; C7 D2 {3 s5 u1 K + j9 F  K: e- f, y
/ x8 ^1 T7 D' b, e# |9 g  z
可以看到,与原曲线图片相比,提取到的数据曲线相似度能达到较高要求。但进一步观察会发现,右图曲线较左图而言,高频分量有一定的减少(即右图曲线更平滑),原因在于数据处理时,对同一个x对应的这组y值进行了均值处理,则在图像上近似反应为均值滤波,从而使得提取到的数据绘制成的曲线的高频分量被滤除。
' t* u/ S9 z5 u7 U* f2 Q3 B+ a( `, [* ~5 K( A- R8 i
最后,将提取到的最终数据,保存起来如下:
; Y4 H  z# p5 b, p3 |! e7 M- m1 W; D0 H0 _6 a/ q% F8 }
3 {# {. m7 ^( c  Z- M$ @- C. d

8 f5 z, S; ~# H) M. P" k/ f1.4进一步的讨论——曲线拟合
# C; p4 {; A! _# M6 ~# d# l
" ^7 s" S4 Y* D2 C7 ]. E通过对图片中曲线的数据提取,可以得到数值上的答案,这会带来进一步的思考,即能否得到这些数据的解析表达式。很容易想到,利用最小二乘法来拟合这些数据,这便涉及到了曲线的拟合。(插值与拟合可以这么理解:对于数据点集,若均落在曲线上,则该曲线为插值曲线,否则为拟合曲线)5 T+ U" Y2 c2 a8 `
6 {& w2 O& W: e& V2 F# Y
对于一些简单的曲线图片(如下),可以考虑用泰勒级数来近似,即多项式拟合。7 h8 R4 n8 m5 _; w6 S& h8 d
1 C, r! b3 B! r1 U. j- H' j

6 a" h& Y3 _+ a8 }8 d% h
. D. e+ z/ l) ?2 Y; ~; A4 D1 A' e* n数据提取并拟合的结果展示如下
1 R& J! R/ l- X  C7 W* i
) l2 j/ |) z! h+ s0 d& x
$ f& R  O1 X  o, [$ T# J
7 w- l$ p4 A7 C2 x7 ^; H同时还能得到拟合多项式的系数
# @% d1 y( S! `  C0 N6 @- k
8 d  v  k* F; R& n, m
' H2 B* q- c; l9 ?- s- U; B/ I
+ Y/ B. N& z3 E& ^. N) R4 S  x2 I从而得到该曲线的多项式(这里采用四阶多项式)表达式为:
. W* g4 `' |; J9 C$ [5 S9 g" B8 ^" B! t- L$ G: i8 t
6 |$ z+ a1 ~, t+ l1 `" C2 d7 u

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

0 Q& j3 G1 q1 F6 v1 z# j* G- ~% u在MATLAB中,提供了cftool工具箱,其提供了拟合与插值的GUI,使用非常方便,直接在命令窗输入cftool即可调用,cftool界面如下所示,其具体使用方法不在此赘述。! @7 ]# v6 m, C  d- g
+ z* i" ~' |8 L) \) _9 ^4 {+ [
8 [' M* X9 i+ g1 H8 K: ]
% o3 r- s$ Q' u, J0 G1 b  q1 y( n1 P% n
2.MATLAB程序
8 H* o" m8 y: `0 B% `MATLAB源代码如下所示,和以往的风格一样,提供了详细的注释8 @& N! V8 b: X

/ X' ^  Z  O! \$ n& N% //提取图片中的曲线数据
' D: [' e# w( K- Nclear,clc,close all* d1 H" v( d; X! m# Z$ Z/ U
%% //图片与曲线间的定标
/ n* o7 i; S- @0 _im=imread('tu1.jpg'); %//读入图片(替换成需要提取曲线的图片)
% H' M$ f: U1 d# N0 Zim=rgb2gray(im); %//灰度变化2 c8 I+ f5 k  L4 n" W4 O
thresh = graythresh(im); %//二值化阈值
5 z( M/ Q! O* ^1 K3 \- A6 L* V4 p2 Fim=im2bw(im,thresh); %//二值化$ {& G% i+ j8 s3 U& T: A% w0 @* T; t
set(0,'defaultfigurecolor','w')
  a6 f1 x. ?- K" f9 i2 aimshow(im) %//显示图片
4 C# W. d  n  O- F. W[y,x]=find(im==0); %//找出图形中的“黑点”的坐标。该坐标是一维数据。
% N, F8 Q7 j5 x& Dy=max(y)-y; %//将屏幕坐标转换为右手系笛卡尔坐标- v2 F# S; |7 J# y( B
y=fliplr(y); %//fliplr()——左右翻转数组& s; z8 w' B, F% o4 a. D
plot(x,y,'r.','Markersize', 2);! s9 U& d, G8 w( p: T4 w! E* `
disp('请在Figrure中先后点击实际坐标框的两个顶点(左上点和右下点),即A、B两点. ');
8 Z5 s* M& r8 _* m[Xx,Yy]=ginput(2); %//Xx,Yy——指实际坐标框的两个顶点
% u+ [) \7 c" }+ M+ Kmin_x=input('最小的x值');  %//输入x轴最小值
1 C' O$ W3 h- ]# T' `/ Dmax_x=input('最大的x值');  %//输入x轴最大值2 `; o, R0 H9 Y% ]9 O" n& |
min_y=input('最小的y值');  %//输入y轴最小值6 u  E5 D; ]( k5 p; m) `% B. N8 G
max_y=input('最大的y值');  %//输入y轴最大值
( @" Z3 V% s; ^5 H  n; sx=(x-Xx(1))*(max_x-min_x)/(Xx(2)-Xx(1))+min_x;5 H0 ?: }8 O! ?& E; N0 n  p
y=(y-Yy(1))*(min_y-max_y)/(Yy(2)-Yy(1))+max_y;; {$ \0 q$ \+ ~( a; n4 [: x& b
plot(x,y,'r.','Markersize', 2);
* u& q3 G2 j$ p5 J( c1 S5 j  ^$ Baxis([min_x,max_x,min_y,max_y])  %//根据输入设置坐标范围
5 N$ ~# i! C, p$ X% [title('由原图片得到的未处理散点图')- p2 s2 B9 Q" [2 X( X' @- Z; R
%% //将散点转换为可用的曲线
) T& s4 Y0 W$ x6 R0 D2 A%//需处理的问题与解决思路
  ?' b% A6 p* }/ f% k%//(1)散点图中可能一个x对应好几个y <---> 保留mean()-std()到mean()+std()之间的y值 并取平均处理
" ?" ]+ D* m; P3 j+ s) j. Q%//(2)曲线的最前端和最后段干扰较大 <---> 去掉曲线整体的前(如5%)和后5%
2 Y3 _3 r* }+ Q& s- N( D' {1 W4 j%//(3)曲线的最顶端和最底段干扰较大 <---> 去掉曲线整体的上10%和下10%; u; y4 m8 m, K8 }. @$ @
& Y$ H4 M8 P, j! P- t$ E" l" C6 X
%//参数预设
" q$ _( l) L/ `2 e( W/ Drate_x=0.08;  %//曲线的最前端和最后段删除比例7 e- `2 I; G3 q7 D4 L8 E' a/ q
rate_y=0.05;  %//曲线的最顶端和最底段删除比例
  c8 n% t8 w9 k0 Y( [9 S- M3 w/ E$ g! P" A
[x_uni,index_x_uni]=unique(x);  %//找出有多少个不同的x坐标
! a( y: m2 E; i6 y( H. L' \1 Q4 Y+ n# D7 ^( ?
x_uni(1:floor(length(x_uni)*rate_x))=[];  %//除去前rate_x(如5%)的x坐标
& J2 E; f) i% |& _- k) qx_uni(floor(length(x_uni)*(1-rate_x)):end)=[];  %//除去后rate_x的x坐标0 F+ K* B7 z. p1 t, \9 o6 S
index_x_uni(1:floor(length(index_x_uni)*rate_x))=[];  %//除去前rate_x的x坐标8 z2 W6 Y4 a; q
index_x_uni(floor(length(index_x_uni)*(1-rate_x)):end)=[];  %//除去后rate_x的x坐标7 K6 s3 Z9 @# n

) r: d- C0 W& D% W  v4 ~8 K[mxu,~]=size(x_uni);' z4 D: n. U$ z. i
[mx,~]=size(x);& W0 @/ \2 p7 c6 \9 L4 {2 U3 S8 m: e2 Q
for ii=1:mxu
# s, z0 O: O8 Z. G7 v9 ]: [    if ii==mxu% d8 k, @; O- a  d& ~0 A) V
        ytemp=y(index_x_uni(ii):mx);
) e7 G! j; w) B    else
! G1 J5 P+ a  p$ u* h        ytemp=y(index_x_uni(ii):index_x_uni(ii+1));
8 i: I% L% H2 M, t% B) w/ s    end' o' n; X$ {; d5 h1 u
    % //删除方差过大的异常点
4 Y  u$ A/ i& o: S& g1 Q8 _1 O* H    threshold1=mean(ytemp)-std(ytemp);  D: g1 z" _' K: t2 @: P
    threshold2=mean(ytemp)+std(ytemp);
, S9 O0 K9 U7 U    ytemp(find(ytemp<threshold1))=[];  %//删除同一个x对应的一段y中的异常点
4 b% W( |& _4 f& o% b: z    ytemp(find(ytemp>threshold2))=[];  ?. k! ]0 r5 a( V/ k5 N# M  V
    % //删除距顶端和底端较近的点  l: Q' j. |. S$ f$ n) X
    thresholdy=(max_y-min_y)*rate_y;  %//y坐标向阈值) y# c1 s! U3 F) g: n
    ytemp(find(ytemp>max_y-thresholdy))=[];  %//删除y轴向距离顶端与底端距离小于rate_y的坐标
( c6 `2 {4 k9 i2 g7 G$ J    ytemp(find(ytemp<min_y+thresholdy))=[];" `, ?* H- Q' u; r3 y& I
    % //剩下的y求均值1 z2 _5 e8 x- x" e) \. }& d
    y_uni(ii)=mean(ytemp);
1 J# a- L" D) B( F7 Tend$ K7 D7 c. E# {- W9 D
% //此时很多x_uni点处对应的y_uni为空,即NAN,要进一步删去这些空点1 e; P+ Y3 I3 C, ]; S1 a3 d5 b# e
x_uni(find(isnan(y_uni)))=[];
+ e4 y4 \: L+ `1 ]6 M- t8 S0 a/ uy_uni(find(isnan(y_uni)))=[];
: d5 C- Z6 n6 Y% n% //画图- o0 w) P6 n* ]  y% \3 g
figure,plot(x_uni,y_uni),title('经处理后得到的扫描曲线')  b& q1 E% ?7 W4 |$ q1 Y
axis([min_x,max_x,min_y,max_y])  %//根据输入设置坐标范围1 }- m+ ^! C- l2 ^; ?
% //将最终提取到的x与y数据保存) L' g; \, R0 T0 A7 J
curve_val(1,:)=x_uni';% O, [8 _% L& b: T" H) ?$ ^5 i* q  g
curve_val(2,:)=y_uni;/ V" x! G4 X/ D4 X" q* u
%% //对提取出的数据进行拟合(按实际情况进行修改)
* f: j0 m0 }/ k# |9 P7 y[p,s]=polyfit(curve_val(1,:),curve_val(2,:),4);  %//多项式拟合(为避免龙格库塔,多项式拟合阶数不宜太高)
5 x" W4 X: m; t( ~[y_fit,DELTA]=polyval(p,x_uni,s);  %//求拟合后多项式在x_uni对应的y_fit值
5 F# e) A4 x+ |figure,plot(x_uni,y_fit),title('拟合后的曲线'): K: E# e9 g4 [$ }
axis([min_x,max_x,min_y,max_y])  %//根据输入设置坐标范围; _+ C% v/ `3 |0 S- h) X1 F

该用户从未签到

2#
发表于 2021-1-20 19:01 | 只看该作者
利用matlab从图片中提取曲线坐标数据
您需要登录后才可以回帖 登录 | 注册

本版积分规则

关闭

推荐内容上一条 /1 下一条

EDA365公众号

关于我们|手机版|EDA365电子论坛网 ( 粤ICP备18020198号-1 )

GMT+8, 2025-11-24 03:18 , Processed in 0.156250 second(s), 26 queries , Gzip On.

深圳市墨知创新科技有限公司

地址:深圳市南山区科技生态园2栋A座805 电话:19926409050

快速回复 返回顶部 返回列表