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

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

[复制链接]

该用户从未签到

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

EDA365欢迎您登录!

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

x

; ?1 C6 [4 u  F2 N$ X1 R目录. v0 S* p/ k5 s
0.引言
* o* ^/ {: N  ]3 ^0 L5 l/ D9 d% A8 a% F1.思路详解与分析3 C! F5 f4 q1 ~# \# e
2.MATLAB程序7 ^6 X5 l; K' b& D+ O
  ~6 V" h/ O' U& M  Q
0.引言
' q- E7 \$ w: ]4 V) J% O  在读文献的时,经常遇到这样的情况:文章里提出的方法好有趣啊,好想拿文中用的数据来试试看看能不能得到相近的结果,可是文中只有根据原始数据绘制的曲线图,没有数据。如下图所示。8 l; U& T2 r0 Z, G6 M, |& w* c" {

0 v9 R8 z$ g% T0 w$ R 6 ?9 P' g  s, G- v/ C* i- A
! Z/ J* V( S5 C, t; M% i/ B
  此时,如果能从文中把这幅图截取下来,输入到一个函数中去,最后能返回从图片中提取到的曲线的坐标数据,岂不美哉。这便是本文的工作。+ U$ h% b% n7 G
4 v3 L* t  m1 D6 {! S
1.思路详解与分析7 c7 y- I4 B! g- ^8 a
1.1准备待提取数据的曲线图片7 V$ |; h4 J: z: j0 K3 o
: |! J: x" k; R% L! U. E
将待提取数据的曲线的图片(如.jpg格式图片),利用 imread 输入到 matlab 中。( Q2 _* C. F6 m, \# T

+ C$ V* G8 X7 H$ M7 c+ S3 h! @1.2曲线图片预处理与数据转换
$ T, H8 c+ V9 }% `; V4 G  m2 v9 R0 t5 d8 i; G/ z
曲线图片预处理步骤的主要工作包含如下:) A6 b% Q! F9 S5 t/ I- ~' b

! M- ?) w# y1 N; e* r, R(1)图像二值化$ L4 c& u& U$ H$ Y4 n$ X
* n* c, a6 H, c8 h% X5 Z
将输入图像进行二值化处理,但分割得到的结果并不全为数据,其中可能还包括坐标轴等干扰点需要去除。
9 f( y; t5 ^. ]( f7 [/ ^
: c& t' `' t: ^" Q(2)获取从图片像素到曲线坐标的定标数据
9 e) \% Y9 ?6 p$ f9 z: @1 }+ `
) s2 V  ~6 O% b' E9 p7 u" m7 n首先,通过ginput()手动从图片中提取到两个像素点,这两个点分别为曲线坐标框的左上角和右下角。5 N" N; y! z2 r' i2 z6 u
# y, \3 z' G5 }6 m
, @5 u8 L, L: g' G1 ?
0 M; v8 \- j$ Q9 _* x  |% _
此时,便获得了曲线在图片上的像素范围
9 t7 A) V. @1 ]+ a6 c
$ A3 k- r, x+ j  a[x_index_min, x_index_max] & [y_index_min, y_index_max]& y" R* k/ [; b

6 m' r: q4 i5 S; I( i! M/ M然后,手动输入实际曲线的数据坐标范围 [x_min, x_max] & [ymin, y_max_]. 如下所示。
% o) `: k3 @9 ?5 t4 c2 R
- R" @0 l1 m$ x# z$ B) m8 U
" x' f2 o1 _0 c+ l1 J6 c7 Q( F4 e' c$ G4 I: g9 p* y
此时,一方面得到了像素坐标,一方面得到了实际坐标。接下来便利用这对数据,将图片中全部的像素坐标转换到实际坐标。8 [+ G/ F2 ?7 T! C& a) R5 d
3 l2 @0 s4 e" v0 O2 O3 g9 \( ~/ M
最终,得到了由图片提取到的数据散点图,如下:2 T& x+ r7 g* L9 p, O9 i

! U4 B. }7 ]# I- \0 y% f- G
3 @/ I1 t7 t7 d( b* ]) m可以看到,此时得到的结果,虽然曲线与所需要的相近,但曲线外的部分,如坐标轴框、坐标轴刻度等也被转换成了数据,还需要进一步的处理。* ]/ [: O4 U0 {  K, m

4 y4 J  D& u7 ?( a# P8 d1.3数据的进一步处理并得到最终曲线# t1 V% v8 w3 C! j7 {+ X
9 E( }1 n( r; j" \
这一步的主要工作是在数据中除去曲线之外的部分(包括坐标轴框、坐标轴刻度等);以及解决一个x数据对应多个y数据的情况。
5 I6 w3 |) G7 g8 s6 \5 t9 Y6 ^( B6 u" h- Z) Y1 F
显然,坐标轴在整幅数据中,均处于边界位置,因此,很容易想到的一种方法是,设定阈值,将距离边界较近的数据直接删除。这里,设定了两个阈值,一个用来限定x方向上的数据,一个用来限定y方向上的数据。比如设定:rate_x=0.08; rate_y=0.05; 意思是阈值设定为曲线最前端8%和最后段8%的数据,曲线最顶端5%和最底端5%的数据。8 s4 e$ k/ V$ k( d! O* u

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

  T" F8 D& w* e" N: P8 ^) [到这里,我们就将数据的处理部分基本完成了,我们将处理后的数据再次绘制成曲线,便可以得到如下
- p% @0 v+ Q- i, j/ @; o( Z2 R5 ^: `/ J/ l5 a
3 E5 g) F" a% d( b0 K6 S
7 p" n7 q' ~) i
对比处理之前的数据,由于限定了范围,因此曲线图片中带来的坐标框等内容转化而来的数据已经被删去。- e5 j5 }- u* ?

. d3 x) Z; ~. V3 B. [将需要提取坐标的曲线图片,和我们提取并处理后的数据绘制的曲线,放在一起对比如下:
* c# W/ {1 A. C/ ?7 @
; f% s1 a) i$ r( G* s( k2 J
) _5 ?# e- ]& y- A) e$ I
+ S$ }( }& G9 E1 u! m- J可以看到,与原曲线图片相比,提取到的数据曲线相似度能达到较高要求。但进一步观察会发现,右图曲线较左图而言,高频分量有一定的减少(即右图曲线更平滑),原因在于数据处理时,对同一个x对应的这组y值进行了均值处理,则在图像上近似反应为均值滤波,从而使得提取到的数据绘制成的曲线的高频分量被滤除。
* r% i  s5 O' t3 r( A9 T
& S% ^' L6 d$ {# R, o' `2 W3 ]1 R最后,将提取到的最终数据,保存起来如下:
/ B4 q; ^0 `1 d2 Q) I9 U& ^
% p5 X1 J  z0 L  O
, m) E" _( V8 D& F, s* t& [8 _' S# C
1.4进一步的讨论——曲线拟合% A7 N! J5 _  k* H
/ r4 S& W* ~. ~* M8 ]
通过对图片中曲线的数据提取,可以得到数值上的答案,这会带来进一步的思考,即能否得到这些数据的解析表达式。很容易想到,利用最小二乘法来拟合这些数据,这便涉及到了曲线的拟合。(插值与拟合可以这么理解:对于数据点集,若均落在曲线上,则该曲线为插值曲线,否则为拟合曲线)
3 U) \) j0 u3 o" ^! m% s
* |; ^- j$ b( {4 O: ?4 E- l% D对于一些简单的曲线图片(如下),可以考虑用泰勒级数来近似,即多项式拟合。. _+ s  L& s2 q' f
, v4 d% v" l  i2 E
1 ?) n6 p' k; m# I( N; ~

* |, r3 i& Z, Y1 u& Z数据提取并拟合的结果展示如下( D1 |4 ?! c7 N& d9 j! k, O2 ~# Y

& w' n9 j* I* e) t& d
" p4 g; ?( s0 T! Z  F
7 X# A) M( F3 r* a' \同时还能得到拟合多项式的系数
- D/ V1 ]: O; Y
& O8 \2 J% `! S6 t7 U$ X2 M/ X
9 Z+ d" N- u" `9 H6 ]3 A
8 B6 T8 u1 j0 N7 M9 B! Y3 \6 N从而得到该曲线的多项式(这里采用四阶多项式)表达式为:- u* P: m& K! D, \6 m" L6 U

4 I/ Z6 h5 y: m$ V9 a8 W
0 b' R- c8 X& \2 w  Y/ F) A
/ ]& i7 ~" Y0 L6 a  S" ~理论上,泰勒级数可以分解任何函数,但实际上,多项式拟合的次数太高,会出现龙格库塔现象,即摆尾现象。因此,多项式拟合的阶数不易过高,一般低于5阶。对于本文最开始的那幅曲线而言,仅5阶的泰勒级数就显得力不从心了,因此,对于这种存在波动剧烈的函数,可以考虑用傅里叶级数进行拟合,或者如果能提供先验知识,可以直接用先验表达式进行拟合。
. M5 H3 V7 I, Z0 V& V" A# c* K
* ]6 Z! w: _3 A在MATLAB中,提供了cftool工具箱,其提供了拟合与插值的GUI,使用非常方便,直接在命令窗输入cftool即可调用,cftool界面如下所示,其具体使用方法不在此赘述。
! H; F" ^; H8 Y2 K: N4 f+ N$ N/ E: y; ?* D9 H5 O3 n; U8 r

: g3 L9 @7 X0 J% Q3 t2 K, b3 {* @9 N; M* M
2.MATLAB程序  _  R  s5 R) W; _
MATLAB源代码如下所示,和以往的风格一样,提供了详细的注释8 E  A  L2 q# Q" n5 z

7 u0 l% k8 e8 B; X) V. ~/ S% //提取图片中的曲线数据
& O: z9 i! b' ^$ \' y# a" pclear,clc,close all* ^+ c3 H8 d( b
%% //图片与曲线间的定标
6 P$ I- A* Y) ^3 e4 C! Vim=imread('tu1.jpg'); %//读入图片(替换成需要提取曲线的图片)2 m0 @& |2 D: a$ j" N7 M% U0 {
im=rgb2gray(im); %//灰度变化$ i: U; x) G- }# p
thresh = graythresh(im); %//二值化阈值
( W4 d$ R0 R; E6 a* Rim=im2bw(im,thresh); %//二值化
- ~, F$ k, V: Y; Zset(0,'defaultfigurecolor','w')) N+ s. h6 @2 n" w6 f$ t; f
imshow(im) %//显示图片8 T0 t3 m+ {. h4 R: W
[y,x]=find(im==0); %//找出图形中的“黑点”的坐标。该坐标是一维数据。
2 \2 u& [+ [+ G5 ?" H4 i5 `y=max(y)-y; %//将屏幕坐标转换为右手系笛卡尔坐标. n: X, w& y4 c% E. j1 F2 u
y=fliplr(y); %//fliplr()——左右翻转数组
6 n6 \: U$ b! c2 ]* Lplot(x,y,'r.','Markersize', 2);' i( Y- Y& m6 k1 `
disp('请在Figrure中先后点击实际坐标框的两个顶点(左上点和右下点),即A、B两点. ');: L* ]9 r  j8 V* t9 u
[Xx,Yy]=ginput(2); %//Xx,Yy——指实际坐标框的两个顶点
, L. y% T7 U5 O1 H# q; J" Nmin_x=input('最小的x值');  %//输入x轴最小值
& @8 m/ M2 S7 {7 G5 jmax_x=input('最大的x值');  %//输入x轴最大值$ z- _/ E, t9 e/ Y/ X# P
min_y=input('最小的y值');  %//输入y轴最小值3 ]1 |& r8 i* _! @2 g, T- C8 a! v
max_y=input('最大的y值');  %//输入y轴最大值( ^, L, u! D. D0 X, }" Z+ |
x=(x-Xx(1))*(max_x-min_x)/(Xx(2)-Xx(1))+min_x;
$ P) _  {9 U- D- K$ `1 Ny=(y-Yy(1))*(min_y-max_y)/(Yy(2)-Yy(1))+max_y;
& l: g( _, I7 e0 {plot(x,y,'r.','Markersize', 2);8 R$ Z- }9 P9 P
axis([min_x,max_x,min_y,max_y])  %//根据输入设置坐标范围
8 x; p! D5 L! U6 F' B1 h5 `title('由原图片得到的未处理散点图')% w& a# B- \: h! B0 B# X6 L
%% //将散点转换为可用的曲线* L2 N6 s, D; n" p4 L# \  Q' h9 J
%//需处理的问题与解决思路
& `+ r& B) j1 K%//(1)散点图中可能一个x对应好几个y <---> 保留mean()-std()到mean()+std()之间的y值 并取平均处理
: {1 \8 G# [1 g) u: M6 ^%//(2)曲线的最前端和最后段干扰较大 <---> 去掉曲线整体的前(如5%)和后5%
$ J3 @3 a! t: `+ O5 H( S! t%//(3)曲线的最顶端和最底段干扰较大 <---> 去掉曲线整体的上10%和下10%9 o; J0 P) w3 G% G* z  k: b5 Y

3 Z) ]( z7 j/ ]%//参数预设
; ]" L: }' n5 y, f  Grate_x=0.08;  %//曲线的最前端和最后段删除比例; `2 j* S7 A' Y& l2 R, a
rate_y=0.05;  %//曲线的最顶端和最底段删除比例
  S' d! Z: ^4 `3 m9 \- I( H0 |- }& I4 B
[x_uni,index_x_uni]=unique(x);  %//找出有多少个不同的x坐标
! L$ v9 B6 R8 Y+ P+ ?) l: f1 V  O. d: C
x_uni(1:floor(length(x_uni)*rate_x))=[];  %//除去前rate_x(如5%)的x坐标
0 ?9 i8 s8 _- B' V; |9 `x_uni(floor(length(x_uni)*(1-rate_x)):end)=[];  %//除去后rate_x的x坐标; I/ k. a8 C( G3 m
index_x_uni(1:floor(length(index_x_uni)*rate_x))=[];  %//除去前rate_x的x坐标' Y: h+ v8 x9 w$ _
index_x_uni(floor(length(index_x_uni)*(1-rate_x)):end)=[];  %//除去后rate_x的x坐标
& V4 e, I' T1 F2 R5 l" Z& a1 L7 ]" k: B) V5 d
[mxu,~]=size(x_uni);/ m/ o# |, B& A- j( C
[mx,~]=size(x);
3 B' t) A( O; v, \% f  ?for ii=1:mxu. x  t+ v/ \) w* `) `
    if ii==mxu
) a/ a- X! S* B& Y        ytemp=y(index_x_uni(ii):mx);
" p! X1 M1 @) ]* f8 C! `    else
1 C8 A" A9 K  J6 _+ m+ k$ W5 @        ytemp=y(index_x_uni(ii):index_x_uni(ii+1));
, ], N5 u" b) c% Y! m    end
# }  W1 f) h. h) K    % //删除方差过大的异常点) B/ n0 s3 C5 ~+ S) J& H
    threshold1=mean(ytemp)-std(ytemp);
/ ~" C; D8 G3 f% E* R$ x5 n& @3 V    threshold2=mean(ytemp)+std(ytemp);" i. n& a6 L( O- P. s
    ytemp(find(ytemp<threshold1))=[];  %//删除同一个x对应的一段y中的异常点
) b, A9 H/ D) G/ \    ytemp(find(ytemp>threshold2))=[];" `7 S: c4 H+ W  _8 _
    % //删除距顶端和底端较近的点
, L4 ?( u, r3 w    thresholdy=(max_y-min_y)*rate_y;  %//y坐标向阈值
6 T" o5 r5 f+ B. O2 Z& N    ytemp(find(ytemp>max_y-thresholdy))=[];  %//删除y轴向距离顶端与底端距离小于rate_y的坐标: L6 J" [( C1 H8 e
    ytemp(find(ytemp<min_y+thresholdy))=[];+ C" s1 i/ `4 ^( }3 ]9 D" O
    % //剩下的y求均值
: H( F6 [. @  U9 e# I4 J9 V    y_uni(ii)=mean(ytemp);
1 J" o2 p* u. X6 w  zend
# Q  [8 Y5 v2 i; ?0 g& _/ E% //此时很多x_uni点处对应的y_uni为空,即NAN,要进一步删去这些空点
- t; c. x! Y; E: K/ S$ Ax_uni(find(isnan(y_uni)))=[];$ i9 m$ z  V! H5 v
y_uni(find(isnan(y_uni)))=[];
. X7 q% V3 I" {# ^% //画图
4 b# D4 N: |- G; f/ k6 n( ]& Q$ Tfigure,plot(x_uni,y_uni),title('经处理后得到的扫描曲线')
6 m0 B4 L8 h6 ?: s5 s" X4 saxis([min_x,max_x,min_y,max_y])  %//根据输入设置坐标范围
$ M2 F3 Y5 y* R) p. {6 U) t0 P% //将最终提取到的x与y数据保存- d& Y! G) P+ E0 Z' i1 r8 |+ p; A
curve_val(1,:)=x_uni';2 G$ _) V$ G9 d& S2 ]6 i
curve_val(2,:)=y_uni;# {1 f  @8 G! w0 A! e
%% //对提取出的数据进行拟合(按实际情况进行修改)6 {1 L  i5 \' L3 k9 K8 U
[p,s]=polyfit(curve_val(1,:),curve_val(2,:),4);  %//多项式拟合(为避免龙格库塔,多项式拟合阶数不宜太高)# G2 r9 b; q) O
[y_fit,DELTA]=polyval(p,x_uni,s);  %//求拟合后多项式在x_uni对应的y_fit值
  Y* y7 q" N7 ~$ Q9 W0 X! vfigure,plot(x_uni,y_fit),title('拟合后的曲线')
6 a4 ]. a9 L0 q/ _& C: ~) l+ Waxis([min_x,max_x,min_y,max_y])  %//根据输入设置坐标范围
  A$ _/ }) Y: a0 j

该用户从未签到

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

本版积分规则

关闭

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

EDA365公众号

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

GMT+8, 2025-8-2 11:23 , Processed in 0.125000 second(s), 26 queries , Gzip On.

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

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

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