|
|
EDA365欢迎您登录!
您需要 登录 才可以下载或查看,没有帐号?注册
x
% U+ i4 U7 u. {目录
, K3 h5 P( y( B; _( {5 Q3 u0.引言
, o, ?, d; i6 i2 Y' R, E1.思路详解与分析
) Y( f& u; V9 U& R2.MATLAB程序
6 ^0 c' k% C# d# Y% A; \
& W1 u' O+ j/ k7 ^0.引言6 j+ z8 f h3 A* Q6 C
在读文献的时,经常遇到这样的情况:文章里提出的方法好有趣啊,好想拿文中用的数据来试试看看能不能得到相近的结果,可是文中只有根据原始数据绘制的曲线图,没有数据。如下图所示。
3 [9 h% j6 b7 X/ \" B8 t* ?
/ B7 m. G; Y+ }2 Y+ Y% d
2 z* h) v( Z, O8 b& }& P
4 l" Y/ k+ K1 C) o# @8 S
此时,如果能从文中把这幅图截取下来,输入到一个函数中去,最后能返回从图片中提取到的曲线的坐标数据,岂不美哉。这便是本文的工作。! e" {0 U: q3 h6 h: j5 C
# W6 }/ h# w: o& W
1.思路详解与分析
' u$ ~4 L8 d; G" X& o" L6 R. X1.1准备待提取数据的曲线图片7 X8 `3 l" V& h' @" z* a1 a
8 ?! B2 V: d! L8 R, X
将待提取数据的曲线的图片(如.jpg格式图片),利用 imread 输入到 matlab 中。
- K& @* m$ s, @6 f5 B0 z7 c5 A; D# q9 _0 d
1.2曲线图片预处理与数据转换/ ]! Q8 X7 x7 K* q
7 F2 Y9 g% s3 K4 G" ~" c
曲线图片预处理步骤的主要工作包含如下:# A- w6 h% C( _9 H& }) S Y
/ I+ a F) L! H3 h, J
(1)图像二值化
3 }( q3 | Y: @( b4 w4 C8 d. L
, t8 c' V3 }4 w, i! W& R, Y0 c将输入图像进行二值化处理,但分割得到的结果并不全为数据,其中可能还包括坐标轴等干扰点需要去除。
% F: Z% [) d$ U Y. g0 }9 M4 A; x
(2)获取从图片像素到曲线坐标的定标数据
) r; k8 F% q: q6 j% y& b' W* m7 Q# m* Q
首先,通过ginput()手动从图片中提取到两个像素点,这两个点分别为曲线坐标框的左上角和右下角。, Y1 g w" g2 d! q
1 O; {1 E& i1 R* q3 @4 f9 g
; X6 g; {4 S# i1 c9 [! A+ Y4 n- K
/ c, n7 P( `. H9 f
此时,便获得了曲线在图片上的像素范围5 t4 j. p* o0 `/ W' k8 M
8 f: e! C- ~% q4 n( ~
[x_index_min, x_index_max] & [y_index_min, y_index_max]. E( p+ I4 T& x& Z3 h7 }
; ^. w" b) _* p) z
然后,手动输入实际曲线的数据坐标范围 [x_min, x_max] & [ymin, y_max_]. 如下所示。
, N. Y- E+ g( d8 `. ^, A
|7 |1 w; E, b- s1 a! V4 [0 r- o& F
/ a @% \7 c6 C1 o% R, }8 q U0 s. r
3 p9 D/ a; W2 f8 P) H0 ?, i2 ^7 `
此时,一方面得到了像素坐标,一方面得到了实际坐标。接下来便利用这对数据,将图片中全部的像素坐标转换到实际坐标。2 e0 g# N2 f( q% O" p9 O# ?$ s
* G8 N0 h: p& M( w+ S2 u7 B/ J
最终,得到了由图片提取到的数据散点图,如下:/ n: e5 M. {: X2 t
. f- }- d4 ^9 b
' n1 E$ C! `1 H2 R# `% t- b
可以看到,此时得到的结果,虽然曲线与所需要的相近,但曲线外的部分,如坐标轴框、坐标轴刻度等也被转换成了数据,还需要进一步的处理。
y1 r) O+ `. @5 B% W% j* d" j4 z5 F: H/ q; G+ X
1.3数据的进一步处理并得到最终曲线2 W1 V, J8 n( d. v
4 [* x1 I$ k. q5 L: X7 [这一步的主要工作是在数据中除去曲线之外的部分(包括坐标轴框、坐标轴刻度等);以及解决一个x数据对应多个y数据的情况。( [. N) v. O! P0 T0 e
8 O8 l; z: }' Y6 U$ M显然,坐标轴在整幅数据中,均处于边界位置,因此,很容易想到的一种方法是,设定阈值,将距离边界较近的数据直接删除。这里,设定了两个阈值,一个用来限定x方向上的数据,一个用来限定y方向上的数据。比如设定:rate_x=0.08; rate_y=0.05; 意思是阈值设定为曲线最前端8%和最后段8%的数据,曲线最顶端5%和最底端5%的数据。% k( v' L7 j9 T8 u) O! o
& g& {/ P: h/ I: J" e! _' S进一步的,对于提取到的数据图,大多数情况一个x会对应若干y,因为数据是由图片转化而来,而图片的分辨率有限,一个实际数据点会用多个像素来表示。解决此问题的中心思想是将同一个x对应的若干个y取均值,但不能直接求均值,因为还有很多y是噪声(如坐标轴线、由图片噪声带入的干扰点等)需要先去除,在第一个问题中,通过限定y的范围,已经在一定程度上除去了部分干扰,在此基础上,我们求取一个x对应的这组y值的均值mean与标准差std,当某些y值位于[mean-std , mean+std]之外,则认为这些y值波动太大,将它们删去。
c& Y- M! E9 e$ Y, b$ |' V1 e/ |1 T- V+ K F0 g
到这里,我们就将数据的处理部分基本完成了,我们将处理后的数据再次绘制成曲线,便可以得到如下
4 j: o1 c; C w4 _* [2 i: r; g+ F+ q9 C
7 |: [8 A- y3 X9 c5 X. I: g
: A" f' b0 u6 x9 F, f' _
对比处理之前的数据,由于限定了范围,因此曲线图片中带来的坐标框等内容转化而来的数据已经被删去。1 [+ f' |' z6 z
. y4 C% d) h8 r- q" M5 i% Y) Q
将需要提取坐标的曲线图片,和我们提取并处理后的数据绘制的曲线,放在一起对比如下:
: {$ G8 R! [9 V4 h! V5 w/ f8 S, u" q; Y: P1 K* o. Y! r' A
7 x5 S2 Q% O# {4 D& Y
9 o7 c+ i7 y& M ~1 V. I
可以看到,与原曲线图片相比,提取到的数据曲线相似度能达到较高要求。但进一步观察会发现,右图曲线较左图而言,高频分量有一定的减少(即右图曲线更平滑),原因在于数据处理时,对同一个x对应的这组y值进行了均值处理,则在图像上近似反应为均值滤波,从而使得提取到的数据绘制成的曲线的高频分量被滤除。
0 d! X6 l7 E5 [! \) `. x5 _! M |2 b2 B& z/ `. P
最后,将提取到的最终数据,保存起来如下:2 r2 H. Q0 s7 k3 A1 f4 p; U
6 H3 e0 b5 I9 {% g- B
! |( q/ [ Y! r/ o- ^, y6 s+ F
: V% E* V+ M& {' n1 D& L. R1.4进一步的讨论——曲线拟合1 ^$ w2 [% z3 ^2 I
^3 p5 R" K3 m% w$ g8 @" W1 Z通过对图片中曲线的数据提取,可以得到数值上的答案,这会带来进一步的思考,即能否得到这些数据的解析表达式。很容易想到,利用最小二乘法来拟合这些数据,这便涉及到了曲线的拟合。(插值与拟合可以这么理解:对于数据点集,若均落在曲线上,则该曲线为插值曲线,否则为拟合曲线)
2 h9 [4 X6 {" V6 o8 Y+ N8 G0 Z+ D
对于一些简单的曲线图片(如下),可以考虑用泰勒级数来近似,即多项式拟合。
; X3 O! U! F9 D7 s
2 p& Y, }* h& R9 ~
, g4 ]8 Z9 X; W7 s {* X) h1 j
( f9 G9 u6 R& Y. W, c
数据提取并拟合的结果展示如下7 p3 S& H& @9 `* @/ j% m, P
# j4 B4 q. @- n# g1 w
0 x5 l/ u6 z- D6 Z1 @! D) [% O t5 v
同时还能得到拟合多项式的系数# n# O. ]6 O+ n g) C/ x. x
& \. p2 z9 @' j6 g* n' {
1 U" M! t+ W; K$ D6 ?/ h9 }( }
+ \' ?; j7 [$ s5 P从而得到该曲线的多项式(这里采用四阶多项式)表达式为:
7 R6 e, C7 j/ d7 I3 ^" W; c% o4 x3 Y
* [0 Y3 a# d' O7 Z6 d H7 t# H9 o9 e! m# {
理论上,泰勒级数可以分解任何函数,但实际上,多项式拟合的次数太高,会出现龙格库塔现象,即摆尾现象。因此,多项式拟合的阶数不易过高,一般低于5阶。对于本文最开始的那幅曲线而言,仅5阶的泰勒级数就显得力不从心了,因此,对于这种存在波动剧烈的函数,可以考虑用傅里叶级数进行拟合,或者如果能提供先验知识,可以直接用先验表达式进行拟合。
% {1 K* b5 ?0 L( ~/ k; W
1 z O! Z) ^- W [8 C在MATLAB中,提供了cftool工具箱,其提供了拟合与插值的GUI,使用非常方便,直接在命令窗输入cftool即可调用,cftool界面如下所示,其具体使用方法不在此赘述。
# ^" @! r+ L, A: F( ~
% s/ R7 e H3 q! m. Z
+ m m2 B8 m( a7 @& O/ y
/ P0 t g# A$ {& O1 N( m a2.MATLAB程序
/ s# A9 ?# y" a* w: ZMATLAB源代码如下所示,和以往的风格一样,提供了详细的注释; X; \4 X. W( R0 [' P0 ]" F
. G }! h/ r; ?' g( {9 }/ A6 m% //提取图片中的曲线数据/ H* `. K% N5 e: J& S( z
clear,clc,close all
/ k7 V# ?6 ?( }$ w%% //图片与曲线间的定标
3 x; T8 [" Z4 z0 \im=imread('tu1.jpg'); %//读入图片(替换成需要提取曲线的图片)* V H* Q2 z1 |+ [; V
im=rgb2gray(im); %//灰度变化& j# e2 |8 P% d, o% e
thresh = graythresh(im); %//二值化阈值1 O' {5 r$ j2 o" J7 J
im=im2bw(im,thresh); %//二值化
9 @) \& D; X1 O; r* { s* C: Sset(0,'defaultfigurecolor','w')1 _0 {0 P! c" P! @
imshow(im) %//显示图片
! F* \! c* d5 n$ l[y,x]=find(im==0); %//找出图形中的“黑点”的坐标。该坐标是一维数据。* q( S, `( v y# {
y=max(y)-y; %//将屏幕坐标转换为右手系笛卡尔坐标: M4 h$ C4 ^' M9 c# @. I
y=fliplr(y); %//fliplr()——左右翻转数组
, {9 E3 K2 X8 g8 ~# Z. m/ P( lplot(x,y,'r.','Markersize', 2);; j j- @2 ~8 o: W1 M
disp('请在Figrure中先后点击实际坐标框的两个顶点(左上点和右下点),即A、B两点. ');! N, d/ a* f: U6 P$ T& a
[Xx,Yy]=ginput(2); %//Xx,Yy——指实际坐标框的两个顶点' _$ E8 Q4 x6 D
min_x=input('最小的x值'); %//输入x轴最小值$ P9 e! u# N9 s' U
max_x=input('最大的x值'); %//输入x轴最大值0 ~: Y9 X4 M2 h" Q
min_y=input('最小的y值'); %//输入y轴最小值; L: j9 g- e- @* P
max_y=input('最大的y值'); %//输入y轴最大值
1 L1 o5 A$ Q' D& N# D' ]' dx=(x-Xx(1))*(max_x-min_x)/(Xx(2)-Xx(1))+min_x;+ A- T1 n% ?# F5 ], n& v d7 a, b
y=(y-Yy(1))*(min_y-max_y)/(Yy(2)-Yy(1))+max_y;
/ X: F7 ^: w7 x- u: l8 eplot(x,y,'r.','Markersize', 2);
: |2 m! J% X9 `/ x! i' M" Faxis([min_x,max_x,min_y,max_y]) %//根据输入设置坐标范围5 R7 I& W# ]3 B1 u3 n1 B; [- o- E. U
title('由原图片得到的未处理散点图')* [2 @6 D2 t* p+ J
%% //将散点转换为可用的曲线& D9 G% b2 F( p0 b4 N0 G
%//需处理的问题与解决思路( f% M: J2 ~7 }8 R7 J2 a
%//(1)散点图中可能一个x对应好几个y <---> 保留mean()-std()到mean()+std()之间的y值 并取平均处理
" W: B; {( [& y8 {%//(2)曲线的最前端和最后段干扰较大 <---> 去掉曲线整体的前(如5%)和后5%# L6 @' a" F# e. A3 X! T8 R$ g
%//(3)曲线的最顶端和最底段干扰较大 <---> 去掉曲线整体的上10%和下10%
" Y/ H: ?4 U v, U4 U, a& l7 w. L$ P% A1 q3 j- M
%//参数预设
0 d& H: N' u: X R2 W( P' Mrate_x=0.08; %//曲线的最前端和最后段删除比例
; b* d- Y3 { o' [" b6 H9 ~rate_y=0.05; %//曲线的最顶端和最底段删除比例" H7 e7 K1 i( [/ Z5 l K
3 X3 g& e- U! p[x_uni,index_x_uni]=unique(x); %//找出有多少个不同的x坐标$ o9 F1 H ]# k, w( v
0 ?6 j' c5 V q+ A' C4 R2 w
x_uni(1:floor(length(x_uni)*rate_x))=[]; %//除去前rate_x(如5%)的x坐标
% _2 x7 ^4 e$ G2 v1 [x_uni(floor(length(x_uni)*(1-rate_x)):end)=[]; %//除去后rate_x的x坐标
+ j% }' p* P! F: P- E$ s$ I# mindex_x_uni(1:floor(length(index_x_uni)*rate_x))=[]; %//除去前rate_x的x坐标
( P: o; z. t; j2 ?0 ]3 xindex_x_uni(floor(length(index_x_uni)*(1-rate_x)):end)=[]; %//除去后rate_x的x坐标7 n" X8 V# R5 W6 n' P$ k
/ U( L! {% X; i4 E* f
[mxu,~]=size(x_uni);
p% |# }2 N5 N5 v, a( K[mx,~]=size(x);
7 W! k0 M+ b/ D2 n- Gfor ii=1:mxu
8 ?- f- {( V# g if ii==mxu4 T" Z3 A& |- j7 ~) y
ytemp=y(index_x_uni(ii):mx);
, U; l, N6 v" d" B7 d V else
! }7 s3 t5 h6 f/ ]' h' N( n: H ytemp=y(index_x_uni(ii):index_x_uni(ii+1));
/ F- x E1 E" M/ w$ s1 w1 a end( F9 y4 i' B8 ]
% //删除方差过大的异常点; l ?9 o) s$ ]4 D* J* F
threshold1=mean(ytemp)-std(ytemp);" t4 @9 ?0 w9 u: W8 L
threshold2=mean(ytemp)+std(ytemp);
8 f1 S% S8 }; U5 Z `) s* W$ ? ytemp(find(ytemp<threshold1))=[]; %//删除同一个x对应的一段y中的异常点
0 }5 T* u! F0 z ytemp(find(ytemp>threshold2))=[];
1 m0 B$ G3 B# D- Y% q- m, n: t, { % //删除距顶端和底端较近的点
4 i3 u% P. |( E# T4 r" W3 a0 c thresholdy=(max_y-min_y)*rate_y; %//y坐标向阈值
( G9 V2 Q5 n4 ?& r( X% N ytemp(find(ytemp>max_y-thresholdy))=[]; %//删除y轴向距离顶端与底端距离小于rate_y的坐标
! r/ J! U4 H( Y ytemp(find(ytemp<min_y+thresholdy))=[];/ n' t( Q8 o% L0 p0 E
% //剩下的y求均值2 g p! W2 U% Q0 Q$ n N
y_uni(ii)=mean(ytemp);
6 f5 {, }6 K9 W/ N9 b5 X3 ?& }2 rend7 @9 @, R3 y) ?. D. A, S2 b2 ^
% //此时很多x_uni点处对应的y_uni为空,即NAN,要进一步删去这些空点
6 r. K5 S, t+ [3 i% l% m: rx_uni(find(isnan(y_uni)))=[];
7 E; J1 Q. D# p0 m5 jy_uni(find(isnan(y_uni)))=[];
4 \+ Y3 _9 b) B% //画图0 u. \5 Q9 A. a" Y* g( n$ Q
figure,plot(x_uni,y_uni),title('经处理后得到的扫描曲线')$ y, i; b' L# j9 W4 W3 k
axis([min_x,max_x,min_y,max_y]) %//根据输入设置坐标范围! ]' D+ K+ W$ n0 W+ U: K; o+ b4 [3 Z
% //将最终提取到的x与y数据保存
) m' N; p- R0 |) Xcurve_val(1,:)=x_uni';1 x p9 `3 ^( G. ?3 I* ^: A
curve_val(2,:)=y_uni;& t: b" E" M5 P+ }4 ^( l& ]& |
%% //对提取出的数据进行拟合(按实际情况进行修改); L6 U7 {' R) b/ a9 d
[p,s]=polyfit(curve_val(1,:),curve_val(2,:),4); %//多项式拟合(为避免龙格库塔,多项式拟合阶数不宜太高)
& u6 Z( y/ G$ N- U1 S8 I `[y_fit,DELTA]=polyval(p,x_uni,s); %//求拟合后多项式在x_uni对应的y_fit值8 p6 U9 |- `! e% b, @! F
figure,plot(x_uni,y_fit),title('拟合后的曲线')
# ] J* d9 c1 C: Z+ N7 s8 uaxis([min_x,max_x,min_y,max_y]) %//根据输入设置坐标范围
4 B1 r4 a' [( Y. ~+ i |
|