|
EDA365欢迎您登录!
您需要 登录 才可以下载或查看,没有帐号?注册
x
! p, o6 j! m) V# |! A5 x
目录
) X) ^. d1 m* H6 p' y0.引言
% I; s' L; U# _$ @1 }. }* ?1.思路详解与分析
2 g" i6 i3 X! Q. ^' W/ f( I2.MATLAB程序
2 D3 G6 I4 L7 c7 t A3 M M. U
+ A E/ F4 H8 `9 Y
0.引言, P6 d; r5 i J8 ]* e W) Z% {
在读文献的时,经常遇到这样的情况:文章里提出的方法好有趣啊,好想拿文中用的数据来试试看看能不能得到相近的结果,可是文中只有根据原始数据绘制的曲线图,没有数据。如下图所示。
1 t8 ^; P3 _4 M% V1 v2 w
# |1 l6 y3 O( \7 o) p% R0 Q7 f4 B% P0 Q3 ]
v+ \% h# A( g8 c% b2 ^
: L& u( C* m1 M p7 f7 o
此时,如果能从文中把这幅图截取下来,输入到一个函数中去,最后能返回从图片中提取到的曲线的坐标数据,岂不美哉。这便是本文的工作。
5 i7 B3 |0 m, ~5 {# K2 U$ V% j9 k& y0 Y: _- E
1.思路详解与分析
0 n* i4 u% K$ P( ]" A- K9 W1.1准备待提取数据的曲线图片
T0 r9 F7 e( I# s! D5 J6 v
; z4 \- P; x) q+ D3 z将待提取数据的曲线的图片(如.jpg格式图片),利用 imread 输入到 matlab 中。
4 t8 Q' A; J7 \# Z# u9 B9 J& x. V: a1 f' \4 F, }
1.2曲线图片预处理与数据转换
- ?* K$ m% Y! `# N6 Z; G7 N2 n* j' z& F, V4 U
曲线图片预处理步骤的主要工作包含如下:
2 H) m# {9 J$ x9 A7 C6 [* y+ ?4 D: p: {1 b4 M4 U8 j' P& r! j
(1)图像二值化, l9 J4 f; b, a3 F3 C
' D! ?% G7 ]* L3 Y7 I
将输入图像进行二值化处理,但分割得到的结果并不全为数据,其中可能还包括坐标轴等干扰点需要去除。. X3 L& I- P: Q# \( h: J
0 [/ b. Y) y$ K+ c! \
(2)获取从图片像素到曲线坐标的定标数据
" n5 R$ {2 b* x9 _
, T4 W5 k- q) C' o C# Z6 x首先,通过ginput()手动从图片中提取到两个像素点,这两个点分别为曲线坐标框的左上角和右下角。
! F1 n* U$ w' h2 R
+ d& z, T0 k. I2 Y+ ?2 F( S! J! O# o
B3 g$ D A; P
4 Y, e" b; \1 G5 ^1 m' f) ]" W
此时,便获得了曲线在图片上的像素范围
6 [4 l3 f+ _9 j
, C" ]6 b+ C9 i, f[x_index_min, x_index_max] & [y_index_min, y_index_max]* y0 E. P \- n6 T
2 i4 {1 a/ B ]: g K然后,手动输入实际曲线的数据坐标范围 [x_min, x_max] & [ymin, y_max_]. 如下所示。+ \- }/ R- N) q- g( S, g
: I5 t- N- y1 U+ o0 V
5 Q1 L( E# h" N; c8 p; s6 h% |/ i M: p w% P9 L; W
此时,一方面得到了像素坐标,一方面得到了实际坐标。接下来便利用这对数据,将图片中全部的像素坐标转换到实际坐标。
5 ~7 I! }) V+ S$ I. ~4 f8 a j
9 l' }+ i7 f/ G( {8 [; K: S最终,得到了由图片提取到的数据散点图,如下:; y+ i8 Y# c) w0 w+ q. L
# E; b1 h5 L, c5 ?& [" Z
, ]6 ~4 B* t% t& v g( c& @9 F
可以看到,此时得到的结果,虽然曲线与所需要的相近,但曲线外的部分,如坐标轴框、坐标轴刻度等也被转换成了数据,还需要进一步的处理。0 f9 t e, j8 Q H) @
2 @) M) H% n: Q, G1.3数据的进一步处理并得到最终曲线( y; i' ]$ M6 {+ c0 z+ R3 _2 J
7 q) g$ {+ M4 E: Z9 W5 Y
这一步的主要工作是在数据中除去曲线之外的部分(包括坐标轴框、坐标轴刻度等);以及解决一个x数据对应多个y数据的情况。% b4 c0 ?& B% Z7 A
' h, w$ u1 z2 r8 Q" u
显然,坐标轴在整幅数据中,均处于边界位置,因此,很容易想到的一种方法是,设定阈值,将距离边界较近的数据直接删除。这里,设定了两个阈值,一个用来限定x方向上的数据,一个用来限定y方向上的数据。比如设定:rate_x=0.08; rate_y=0.05; 意思是阈值设定为曲线最前端8%和最后段8%的数据,曲线最顶端5%和最底端5%的数据。$ y- t2 @) A1 I
9 `7 S. A& o' A6 S$ |2 A5 c进一步的,对于提取到的数据图,大多数情况一个x会对应若干y,因为数据是由图片转化而来,而图片的分辨率有限,一个实际数据点会用多个像素来表示。解决此问题的中心思想是将同一个x对应的若干个y取均值,但不能直接求均值,因为还有很多y是噪声(如坐标轴线、由图片噪声带入的干扰点等)需要先去除,在第一个问题中,通过限定y的范围,已经在一定程度上除去了部分干扰,在此基础上,我们求取一个x对应的这组y值的均值mean与标准差std,当某些y值位于[mean-std , mean+std]之外,则认为这些y值波动太大,将它们删去。
' @( {7 A3 d- R8 R( _! O& ]6 B0 a+ c, Z
到这里,我们就将数据的处理部分基本完成了,我们将处理后的数据再次绘制成曲线,便可以得到如下1 v$ K/ [) t; b& B, Z& v
6 ], ^: V2 F/ [& f. _- r: B
/ q( {3 L4 B3 ^# |
2 K" ]# o* S$ o1 k
对比处理之前的数据,由于限定了范围,因此曲线图片中带来的坐标框等内容转化而来的数据已经被删去。
6 {) a! s& j( P: r* X# ^. d. }+ C( G) O* K& U% Y% f
将需要提取坐标的曲线图片,和我们提取并处理后的数据绘制的曲线,放在一起对比如下:
/ M% D$ {3 O s. T
( f: U+ |: s% \. b3 {: q' m
9 v7 |+ v7 _: C1 _
# D6 S! ?4 n3 b. {* }8 X" ^可以看到,与原曲线图片相比,提取到的数据曲线相似度能达到较高要求。但进一步观察会发现,右图曲线较左图而言,高频分量有一定的减少(即右图曲线更平滑),原因在于数据处理时,对同一个x对应的这组y值进行了均值处理,则在图像上近似反应为均值滤波,从而使得提取到的数据绘制成的曲线的高频分量被滤除。8 E0 G4 D" b2 A' g1 a
8 v2 S* U0 }6 r* {/ n
最后,将提取到的最终数据,保存起来如下:6 F! V$ n& K! \! ~% j* m* j8 \
2 S/ m2 r& N1 T% u5 t4 K; @# R" V) R
6 o4 f* K* W$ m# p& b% f5 ~( m' q/ R/ l' {3 @7 h
1.4进一步的讨论——曲线拟合/ n5 Y9 _0 X+ X6 @ D; z
! M$ x0 o2 _7 {9 r- _
通过对图片中曲线的数据提取,可以得到数值上的答案,这会带来进一步的思考,即能否得到这些数据的解析表达式。很容易想到,利用最小二乘法来拟合这些数据,这便涉及到了曲线的拟合。(插值与拟合可以这么理解:对于数据点集,若均落在曲线上,则该曲线为插值曲线,否则为拟合曲线). f6 V) k( G+ D x
, U8 `( r5 x5 w* N! A对于一些简单的曲线图片(如下),可以考虑用泰勒级数来近似,即多项式拟合。* o& [& {3 z. f
$ x: Z" U/ J: P# z1 r5 q
) F2 b- m3 `4 ^( V
, h2 v! S% t; m" G( I数据提取并拟合的结果展示如下' e$ R' u2 _6 g% l8 U ^' X
H$ p; s1 r! K- K0 Y# \: j
4 D( x: e1 F7 e/ k+ F. H5 q7 _7 X7 P* d
同时还能得到拟合多项式的系数
. w# e4 N$ N7 t2 n, a) ]% U, c. [: H7 E$ P
8 k8 s3 \! @* S* X- {
/ e( o' Z, H1 J+ n6 C7 t0 B从而得到该曲线的多项式(这里采用四阶多项式)表达式为:. }$ e- X o8 h/ G7 V2 K5 h
) c" l/ w. X4 L6 B
" b6 u5 F; ~+ S! F6 }6 l7 q) ^+ n( m( y( G4 _) y+ k8 ?. V2 C
理论上,泰勒级数可以分解任何函数,但实际上,多项式拟合的次数太高,会出现龙格库塔现象,即摆尾现象。因此,多项式拟合的阶数不易过高,一般低于5阶。对于本文最开始的那幅曲线而言,仅5阶的泰勒级数就显得力不从心了,因此,对于这种存在波动剧烈的函数,可以考虑用傅里叶级数进行拟合,或者如果能提供先验知识,可以直接用先验表达式进行拟合。
% B3 W% y8 ^" [- w+ q: i1 b
7 F6 Q8 ^. ^3 D$ x在MATLAB中,提供了cftool工具箱,其提供了拟合与插值的GUI,使用非常方便,直接在命令窗输入cftool即可调用,cftool界面如下所示,其具体使用方法不在此赘述。
: k) v$ @1 n7 H0 ^' `- n* [2 o! _
% b! O b0 Q6 T6 i, j6 w* m
! w1 z6 X) r2 T3 _
2.MATLAB程序# g7 @7 j' P6 T
MATLAB源代码如下所示,和以往的风格一样,提供了详细的注释4 t* F; K& `1 e/ ]- X
) }( L' y; j' D# Y4 h# x
% //提取图片中的曲线数据
$ o5 K' H# K4 _ Z+ M/ Z, Oclear,clc,close all
& R( m- b9 ?7 z% J, X%% //图片与曲线间的定标! q( n" v2 Z" O R0 {7 N0 N' u
im=imread('tu1.jpg'); %//读入图片(替换成需要提取曲线的图片)) ]4 p$ u4 P9 S8 o$ D7 D
im=rgb2gray(im); %//灰度变化3 V' ^2 f7 n' F0 \( ?2 J
thresh = graythresh(im); %//二值化阈值
0 J6 m$ x+ c: U7 l7 s3 Tim=im2bw(im,thresh); %//二值化
% C1 s5 y8 n# j) Jset(0,'defaultfigurecolor','w')- K* y9 R+ M( K P3 u
imshow(im) %//显示图片
. M/ k$ P# H+ T# ~& z[y,x]=find(im==0); %//找出图形中的“黑点”的坐标。该坐标是一维数据。/ J( x2 ]. N( U$ A
y=max(y)-y; %//将屏幕坐标转换为右手系笛卡尔坐标
: v$ {. x& `( G( P9 Y" n; Vy=fliplr(y); %//fliplr()——左右翻转数组
% ]$ l, k: R! r3 oplot(x,y,'r.','Markersize', 2);1 u4 X% D. {; i" n4 h
disp('请在Figrure中先后点击实际坐标框的两个顶点(左上点和右下点),即A、B两点. ');3 s9 y4 W8 H4 Z; @! {0 l
[Xx,Yy]=ginput(2); %//Xx,Yy——指实际坐标框的两个顶点
$ E: ]+ r$ [. ~min_x=input('最小的x值'); %//输入x轴最小值) d5 f' g$ l) Z
max_x=input('最大的x值'); %//输入x轴最大值: q H9 [( V! B' |& Q
min_y=input('最小的y值'); %//输入y轴最小值8 g3 @" Z, s& r6 o4 w
max_y=input('最大的y值'); %//输入y轴最大值0 f- A, n0 t6 F1 Z1 c# V
x=(x-Xx(1))*(max_x-min_x)/(Xx(2)-Xx(1))+min_x;% b7 `0 N M d7 N3 L; J" [
y=(y-Yy(1))*(min_y-max_y)/(Yy(2)-Yy(1))+max_y;
5 L0 m' R4 G( Q9 Fplot(x,y,'r.','Markersize', 2);$ W/ e# f* H& W; Z) E9 `, s
axis([min_x,max_x,min_y,max_y]) %//根据输入设置坐标范围) B9 `6 I9 T) h% n- L+ J7 e, f
title('由原图片得到的未处理散点图')$ @/ ]# H7 B! }' l& n P& x4 x+ k
%% //将散点转换为可用的曲线$ Q- A) O% Q/ \( R
%//需处理的问题与解决思路
/ L# K9 Z( K) c0 R8 o; N) ^0 x9 b%//(1)散点图中可能一个x对应好几个y <---> 保留mean()-std()到mean()+std()之间的y值 并取平均处理
) b2 Z" T( ~5 o t%//(2)曲线的最前端和最后段干扰较大 <---> 去掉曲线整体的前(如5%)和后5%
* ~: \$ K a+ ~ ~: c%//(3)曲线的最顶端和最底段干扰较大 <---> 去掉曲线整体的上10%和下10%/ u$ D( }6 [: X. y
; g6 U- c" v8 p, A6 R5 A
%//参数预设
3 }; Q$ ~1 T0 P/ E9 Lrate_x=0.08; %//曲线的最前端和最后段删除比例
, D( i. X- J# h0 [rate_y=0.05; %//曲线的最顶端和最底段删除比例1 ^2 M# i0 q$ a4 ` m
5 ]- _% T/ i4 \3 B/ d1 t: g% I6 M D
[x_uni,index_x_uni]=unique(x); %//找出有多少个不同的x坐标
4 v% o5 L+ K; |4 [
- P4 Z0 k* n" K( K5 }x_uni(1:floor(length(x_uni)*rate_x))=[]; %//除去前rate_x(如5%)的x坐标: r$ D2 S; ]! c% {) D
x_uni(floor(length(x_uni)*(1-rate_x)):end)=[]; %//除去后rate_x的x坐标
/ F/ Z& G1 A, L. g8 E1 ?: K; ^. uindex_x_uni(1:floor(length(index_x_uni)*rate_x))=[]; %//除去前rate_x的x坐标% B3 x9 I( d6 A, I
index_x_uni(floor(length(index_x_uni)*(1-rate_x)):end)=[]; %//除去后rate_x的x坐标
- L3 z1 o' V& ]0 J3 i7 w4 k; H6 c( }5 l% w7 V2 }5 w8 z4 h" R
[mxu,~]=size(x_uni);
' B6 M4 k1 L1 w' m" T( t[mx,~]=size(x);$ v4 z# B( z, b* ]; S( n
for ii=1:mxu
, m$ s% H t1 r if ii==mxu6 W# f; E* t2 q; N: n- |
ytemp=y(index_x_uni(ii):mx);" X* ~( _- t$ L: L2 a
else
/ X& ^6 `4 u5 g2 A ytemp=y(index_x_uni(ii):index_x_uni(ii+1));
* g/ x1 \. t* c; Y end
7 m+ ~" \- m: n: j % //删除方差过大的异常点
$ s& L0 B/ a+ w; g- T8 Q threshold1=mean(ytemp)-std(ytemp);. ~7 t) O" ^4 R
threshold2=mean(ytemp)+std(ytemp);: O9 {1 j9 N3 W' A3 y
ytemp(find(ytemp<threshold1))=[]; %//删除同一个x对应的一段y中的异常点/ n4 a- v* q( i1 R7 P3 U& \
ytemp(find(ytemp>threshold2))=[];
, ^- D7 h. g/ `/ g% V % //删除距顶端和底端较近的点; G' y! y. B8 e- u' `% q# a& c
thresholdy=(max_y-min_y)*rate_y; %//y坐标向阈值
4 Z6 K# c) r( t; o7 M9 n ytemp(find(ytemp>max_y-thresholdy))=[]; %//删除y轴向距离顶端与底端距离小于rate_y的坐标
7 @1 X4 j# T1 V) ^! } ytemp(find(ytemp<min_y+thresholdy))=[];) z# _: ]1 g/ h" F
% //剩下的y求均值! d7 h3 S& n! u& T2 _2 W
y_uni(ii)=mean(ytemp);; ]6 i; m) C0 f% I. A1 C: a6 c, F6 T
end
2 o5 ?# Z6 S4 K& G; c% //此时很多x_uni点处对应的y_uni为空,即NAN,要进一步删去这些空点
, o C# f* \ c5 xx_uni(find(isnan(y_uni)))=[];
: o, n+ h+ }$ K A9 c* v0 e8 f2 jy_uni(find(isnan(y_uni)))=[];
" v+ b# F- O7 ^1 ]8 s$ w% //画图, x9 N5 ^0 c0 `8 s) g
figure,plot(x_uni,y_uni),title('经处理后得到的扫描曲线')
7 e5 S7 I, A$ }# Z5 ?- kaxis([min_x,max_x,min_y,max_y]) %//根据输入设置坐标范围
/ w( c* D. ^/ n9 u% //将最终提取到的x与y数据保存) L5 y' h. T9 f! t
curve_val(1,:)=x_uni';
^' Y9 X- t8 b, _: g. N2 Fcurve_val(2,:)=y_uni;7 p2 L6 ?& _0 T8 Y u' X- k
%% //对提取出的数据进行拟合(按实际情况进行修改)1 O0 ?' c( f8 e8 I, F
[p,s]=polyfit(curve_val(1,:),curve_val(2,:),4); %//多项式拟合(为避免龙格库塔,多项式拟合阶数不宜太高)
: b$ E8 v* w. j[y_fit,DELTA]=polyval(p,x_uni,s); %//求拟合后多项式在x_uni对应的y_fit值0 K4 f* }. X8 \' A0 E$ p+ B
figure,plot(x_uni,y_fit),title('拟合后的曲线')
: ?# }+ v6 B% \! p+ U$ vaxis([min_x,max_x,min_y,max_y]) %//根据输入设置坐标范围
" y2 \- m$ a1 R* Y+ F) w |
|