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

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

[复制链接]

该用户从未签到

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

EDA365欢迎您登录!

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

x

. G/ K* g" U# r目录: [5 p( d1 x+ d! I  W3 }. S/ \
0.引言# ?7 {' }, a. o. O/ U1 g- [! b
1.思路详解与分析
# n7 ]1 H3 ~4 d2.MATLAB程序5 @! Y- P' ]  X' }, |5 s; w2 W; W

; _) f' F  b" V0.引言( F4 z( Z: X4 D5 q1 Q( j$ U
  在读文献的时,经常遇到这样的情况:文章里提出的方法好有趣啊,好想拿文中用的数据来试试看看能不能得到相近的结果,可是文中只有根据原始数据绘制的曲线图,没有数据。如下图所示。1 S: O1 B' I! u% D) u+ t) g
9 H# U  N; A; y1 h7 G
, k% x. t8 F' N/ i

1 j, p$ y+ l, L( N& ?4 G7 v( Q  此时,如果能从文中把这幅图截取下来,输入到一个函数中去,最后能返回从图片中提取到的曲线的坐标数据,岂不美哉。这便是本文的工作。
# p$ P* N% \* @; q! \5 R2 U! U: O1 u. Q; K) M
1.思路详解与分析
4 R& z& v4 G  E) z1.1准备待提取数据的曲线图片
9 C& L4 a. V0 r3 G3 f9 z
$ O2 l( F1 f% x* l/ ]. B. S将待提取数据的曲线的图片(如.jpg格式图片),利用 imread 输入到 matlab 中。  U; \. O7 o$ u& t, D7 O" i
2 n' t' u5 g: N& ?
1.2曲线图片预处理与数据转换" D( x) x: \8 I7 \- x

  f+ M5 v3 U  H/ U/ ^. z4 f7 S' x曲线图片预处理步骤的主要工作包含如下:
( G# K6 i7 x. T6 ~* \, B5 z5 F+ E1 x- a6 V. ]5 N+ F: x9 E
(1)图像二值化
8 m: w( B6 G3 P& S& l# e) x* U  l3 M; m/ L# V
将输入图像进行二值化处理,但分割得到的结果并不全为数据,其中可能还包括坐标轴等干扰点需要去除。
8 T8 C! N5 U* S- |  b
8 i8 p5 K& r% b) U. E8 A(2)获取从图片像素到曲线坐标的定标数据' _8 a8 X7 a9 w& T2 i

1 M" C  c2 u% N首先,通过ginput()手动从图片中提取到两个像素点,这两个点分别为曲线坐标框的左上角和右下角。
" I; P% f* `+ A2 A, \1 j1 G( y" P; m  H1 z( D' o. U6 E7 D5 a
8 z+ _: ?8 q2 g
  I2 A! E* v0 T5 U3 W
此时,便获得了曲线在图片上的像素范围6 W% W" d0 s3 B' [$ A1 w$ G+ Q/ Y
4 p8 W; e1 i& T8 J# U
[x_index_min, x_index_max] & [y_index_min, y_index_max]
1 x$ B8 E) A5 E
4 F! Q: J  K1 M: S+ V" I* F+ C然后,手动输入实际曲线的数据坐标范围 [x_min, x_max] & [ymin, y_max_]. 如下所示。
4 y: L, x( t* D7 ^3 m6 p7 g( o9 p7 }
/ Z1 H. Z0 Q, }; g' E$ T
$ g( W* J) l0 [/ ^% G+ W& a$ _% |# d) l2 V! H
此时,一方面得到了像素坐标,一方面得到了实际坐标。接下来便利用这对数据,将图片中全部的像素坐标转换到实际坐标。
4 f0 n) h2 j; p" s6 N! j* U8 \$ N: _0 K2 J' U4 [% J' f8 Q
最终,得到了由图片提取到的数据散点图,如下:
- i. d) P. L9 P! n" X8 K$ W( v! D7 M% q8 a4 c4 g0 X, O; G3 @
" v6 v, M* H4 o. h
可以看到,此时得到的结果,虽然曲线与所需要的相近,但曲线外的部分,如坐标轴框、坐标轴刻度等也被转换成了数据,还需要进一步的处理。
* u$ v1 U& T7 r  L  r% R; b8 h6 K. D3 Y
1.3数据的进一步处理并得到最终曲线
2 ^5 D+ Q: f- B
& L  Q% _# x2 r8 R3 `. a5 c这一步的主要工作是在数据中除去曲线之外的部分(包括坐标轴框、坐标轴刻度等);以及解决一个x数据对应多个y数据的情况。
5 A  T+ H! K  c( a; @2 ?/ V( X& F( ]: K2 f
显然,坐标轴在整幅数据中,均处于边界位置,因此,很容易想到的一种方法是,设定阈值,将距离边界较近的数据直接删除。这里,设定了两个阈值,一个用来限定x方向上的数据,一个用来限定y方向上的数据。比如设定:rate_x=0.08; rate_y=0.05; 意思是阈值设定为曲线最前端8%和最后段8%的数据,曲线最顶端5%和最底端5%的数据。: \2 a" P% F, i, n! X
# q! o2 X* g- B. ^! O8 `
进一步的,对于提取到的数据图,大多数情况一个x会对应若干y,因为数据是由图片转化而来,而图片的分辨率有限,一个实际数据点会用多个像素来表示。解决此问题的中心思想是将同一个x对应的若干个y取均值,但不能直接求均值,因为还有很多y是噪声(如坐标轴线、由图片噪声带入的干扰点等)需要先去除,在第一个问题中,通过限定y的范围,已经在一定程度上除去了部分干扰,在此基础上,我们求取一个x对应的这组y值的均值mean与标准差std,当某些y值位于[mean-std , mean+std]之外,则认为这些y值波动太大,将它们删去。" E9 ]: J7 r, A' O3 t/ e' L; x

7 @3 c* v/ ]( z6 ?9 Z4 a到这里,我们就将数据的处理部分基本完成了,我们将处理后的数据再次绘制成曲线,便可以得到如下; L  }4 f# Z* N  ~
" t% c) ]( r# O7 i. R2 @
8 {$ V+ e2 _9 N( l5 O. b
0 l4 M. C+ O$ M; c* b: F% O
对比处理之前的数据,由于限定了范围,因此曲线图片中带来的坐标框等内容转化而来的数据已经被删去。2 b( g, x, v7 A, n4 P" p# j
5 P' h3 Y5 e% p  w4 f( x$ _
将需要提取坐标的曲线图片,和我们提取并处理后的数据绘制的曲线,放在一起对比如下:6 i- e: A) w5 x/ |. ]1 q
* `3 A' t) F  _+ e

( S  T6 ^% X# Q$ k' j
1 u" I  k6 S1 X) `/ m2 K可以看到,与原曲线图片相比,提取到的数据曲线相似度能达到较高要求。但进一步观察会发现,右图曲线较左图而言,高频分量有一定的减少(即右图曲线更平滑),原因在于数据处理时,对同一个x对应的这组y值进行了均值处理,则在图像上近似反应为均值滤波,从而使得提取到的数据绘制成的曲线的高频分量被滤除。6 P# `  F% [4 m3 @
9 i9 t& `$ K! u
最后,将提取到的最终数据,保存起来如下:
5 ?8 f: R6 S( Y9 V. \8 {7 s7 ], J1 a0 a: m4 Z: F: K3 w0 s) r" h* r2 l
& Y0 g# a# I& g& U

% v  E# r4 _; k" l- A! l  b1.4进一步的讨论——曲线拟合
) H4 F1 S7 ?7 L5 x7 a( S) ]+ F
% l1 I+ y  ?8 Y通过对图片中曲线的数据提取,可以得到数值上的答案,这会带来进一步的思考,即能否得到这些数据的解析表达式。很容易想到,利用最小二乘法来拟合这些数据,这便涉及到了曲线的拟合。(插值与拟合可以这么理解:对于数据点集,若均落在曲线上,则该曲线为插值曲线,否则为拟合曲线)* Q1 W$ }' V; j, v; A- r

- {* P4 Q( f, P8 W  l' _对于一些简单的曲线图片(如下),可以考虑用泰勒级数来近似,即多项式拟合。, r# z* Y5 w% z  v) V

4 l* \. W6 P2 p& p  ]# i
, c8 C% [; y' Z- L! W0 u( I6 B7 j4 D9 L8 f3 `" g3 I% ?+ o
数据提取并拟合的结果展示如下/ o2 W* D+ _6 K2 l

; @+ x7 {9 M% n " t" V' D( l! j& ]9 \8 f7 ^) u
& W' N$ \8 F$ Z; p6 ?9 S) v/ y
同时还能得到拟合多项式的系数
$ J. U! r3 n6 L) V
; N# x* `0 l& R' x
* }" C) ]* Q, P& P$ c9 ]4 u1 [# ~6 B0 m& _7 \
从而得到该曲线的多项式(这里采用四阶多项式)表达式为:9 u& |4 j/ j) b+ Q2 Y' I( s5 K

+ ?* `" w& }9 `8 F2 w
0 p. c' V( [1 l9 j! W/ o" L9 q) j; o9 G
' K8 q: l3 K" R/ O理论上,泰勒级数可以分解任何函数,但实际上,多项式拟合的次数太高,会出现龙格库塔现象,即摆尾现象。因此,多项式拟合的阶数不易过高,一般低于5阶。对于本文最开始的那幅曲线而言,仅5阶的泰勒级数就显得力不从心了,因此,对于这种存在波动剧烈的函数,可以考虑用傅里叶级数进行拟合,或者如果能提供先验知识,可以直接用先验表达式进行拟合。
# H6 O" i- Q  k, B# k& }& d" J: H) L+ P/ B7 M1 ?4 _! y
在MATLAB中,提供了cftool工具箱,其提供了拟合与插值的GUI,使用非常方便,直接在命令窗输入cftool即可调用,cftool界面如下所示,其具体使用方法不在此赘述。
: X0 O, U, a; G7 P8 C8 x, c
( ^2 O$ @( z& q- X7 L! g
" i# l7 g+ ~; r# M8 p. q  B
% n; B; U2 @' A3 r: P! O7 i% U2.MATLAB程序. s0 r; o7 c' D5 F$ U
MATLAB源代码如下所示,和以往的风格一样,提供了详细的注释
, Y2 C9 j* g  j, B" @2 R9 B/ a& [# a# W  I0 @1 D* z" z
% //提取图片中的曲线数据
7 h: B7 n$ s1 O) U9 ~$ cclear,clc,close all
7 x2 u  t9 Y6 `: w/ X%% //图片与曲线间的定标
4 i- [# G: `4 ^5 U! P- gim=imread('tu1.jpg'); %//读入图片(替换成需要提取曲线的图片)
3 ~) e7 O. _3 z/ O4 L/ pim=rgb2gray(im); %//灰度变化: w: y) q9 o0 ~: J. k0 Z
thresh = graythresh(im); %//二值化阈值) R4 n& s5 Z' T) N+ `/ x
im=im2bw(im,thresh); %//二值化; m2 ~1 j1 K& ?2 b$ A' t% A/ u2 x
set(0,'defaultfigurecolor','w')
% r  }' w+ b( v7 Q0 d. kimshow(im) %//显示图片
" t+ p6 H2 U7 ~6 l! K' N[y,x]=find(im==0); %//找出图形中的“黑点”的坐标。该坐标是一维数据。1 m4 x6 z4 @) [
y=max(y)-y; %//将屏幕坐标转换为右手系笛卡尔坐标
+ X3 m9 k6 i0 g- E" s2 Y# Gy=fliplr(y); %//fliplr()——左右翻转数组# Y6 q/ e- j6 H" `  X! T
plot(x,y,'r.','Markersize', 2);& t' j; J: ~" i6 Y/ l$ ^& \+ n
disp('请在Figrure中先后点击实际坐标框的两个顶点(左上点和右下点),即A、B两点. ');
- U4 b* D1 g' F( C' C. l[Xx,Yy]=ginput(2); %//Xx,Yy——指实际坐标框的两个顶点
& [8 N  k& M" L9 n% ^2 _1 _min_x=input('最小的x值');  %//输入x轴最小值
3 w) k, P7 ^" p2 gmax_x=input('最大的x值');  %//输入x轴最大值2 u+ u' H9 e, U) v
min_y=input('最小的y值');  %//输入y轴最小值
" _& T) {8 T+ J: g; Amax_y=input('最大的y值');  %//输入y轴最大值
! Y( r* a, f0 ix=(x-Xx(1))*(max_x-min_x)/(Xx(2)-Xx(1))+min_x;
, G" K4 i. F8 B9 v# v: gy=(y-Yy(1))*(min_y-max_y)/(Yy(2)-Yy(1))+max_y;
/ n5 C  K/ @% n4 Rplot(x,y,'r.','Markersize', 2);. J" \5 O8 l' b' M
axis([min_x,max_x,min_y,max_y])  %//根据输入设置坐标范围
5 Z7 P2 O, Q- G; X" m1 b) r3 htitle('由原图片得到的未处理散点图')
7 x4 M8 O- L2 a' Z%% //将散点转换为可用的曲线
/ O, V  t7 H2 ~" ~( V%//需处理的问题与解决思路' a" K, V: K' s( S8 {; d
%//(1)散点图中可能一个x对应好几个y <---> 保留mean()-std()到mean()+std()之间的y值 并取平均处理
; ^& _8 j' Y8 I' F* o%//(2)曲线的最前端和最后段干扰较大 <---> 去掉曲线整体的前(如5%)和后5%
( e3 b6 X( Y  |: G9 T8 Q9 z2 p  y%//(3)曲线的最顶端和最底段干扰较大 <---> 去掉曲线整体的上10%和下10%% Q" ?: g$ R2 C2 m
5 t4 ~- }' z, `4 u4 \) k: C) A- P
%//参数预设
! b1 W$ Z$ h0 m& {) B% Grate_x=0.08;  %//曲线的最前端和最后段删除比例
0 M2 D8 P" ]+ S) Vrate_y=0.05;  %//曲线的最顶端和最底段删除比例8 s& {! |( g: m3 \
) V3 Z+ d* r  {' Z' Q8 D: [0 }1 ]
[x_uni,index_x_uni]=unique(x);  %//找出有多少个不同的x坐标% E1 @( e  ?0 c; h( \! `& `

4 z5 i" [2 I* o1 Wx_uni(1:floor(length(x_uni)*rate_x))=[];  %//除去前rate_x(如5%)的x坐标+ Y* j" j+ {' T7 k2 s6 Z
x_uni(floor(length(x_uni)*(1-rate_x)):end)=[];  %//除去后rate_x的x坐标
0 y6 L; q. K* r, f8 Uindex_x_uni(1:floor(length(index_x_uni)*rate_x))=[];  %//除去前rate_x的x坐标
# w3 {# j1 }( [5 Z- _index_x_uni(floor(length(index_x_uni)*(1-rate_x)):end)=[];  %//除去后rate_x的x坐标
' P5 C7 r8 s; G, O$ k. ~, a4 E+ z$ Z
[mxu,~]=size(x_uni);2 M2 O+ h3 D& w
[mx,~]=size(x);9 j- K, \1 X2 D# x
for ii=1:mxu
7 b1 \+ M! {- ^7 Y- t    if ii==mxu5 w) e7 _# ?/ O/ r, t
        ytemp=y(index_x_uni(ii):mx);
  N/ D1 x) o  N1 u1 Q    else
4 N' z) E8 _) x$ [* ^% u        ytemp=y(index_x_uni(ii):index_x_uni(ii+1));
# F  j3 @6 ~. U    end
" t; ?% G1 {( H3 Z, B* U+ \    % //删除方差过大的异常点
  o( t4 p+ }) o4 c    threshold1=mean(ytemp)-std(ytemp);4 w8 E! X' H5 l  b: S: a% Q' _
    threshold2=mean(ytemp)+std(ytemp);
. }6 i4 `7 s7 W3 ?' a" ^, k' _    ytemp(find(ytemp<threshold1))=[];  %//删除同一个x对应的一段y中的异常点7 x4 q# a$ ?4 u, O; O
    ytemp(find(ytemp>threshold2))=[];# Z! w# e2 U$ d4 s$ M* E; l
    % //删除距顶端和底端较近的点
+ R3 u6 Y# e, D' v% _' I3 o    thresholdy=(max_y-min_y)*rate_y;  %//y坐标向阈值
- M! ^8 Q1 v0 n8 O- C# B. R% H    ytemp(find(ytemp>max_y-thresholdy))=[];  %//删除y轴向距离顶端与底端距离小于rate_y的坐标  V1 ~" A: ?3 J& i* c
    ytemp(find(ytemp<min_y+thresholdy))=[];" y9 W+ j% f/ P& D' f2 T
    % //剩下的y求均值
& ^4 W9 r; i" u    y_uni(ii)=mean(ytemp);
! M: o2 ~; A$ {) d* bend
8 p. P& C  d, v+ t/ R) w! s2 m% //此时很多x_uni点处对应的y_uni为空,即NAN,要进一步删去这些空点- H' w; b" y$ m0 M
x_uni(find(isnan(y_uni)))=[];2 i: I% v) e" t/ r0 c+ T$ x
y_uni(find(isnan(y_uni)))=[];
9 B2 k- e4 p1 S7 h' f7 [% l2 J) }% //画图! ^+ ]2 l0 G4 P8 m9 y
figure,plot(x_uni,y_uni),title('经处理后得到的扫描曲线')
$ o. s! L/ c9 l1 D) qaxis([min_x,max_x,min_y,max_y])  %//根据输入设置坐标范围
% `3 N. e. Z- @% p1 d! z% //将最终提取到的x与y数据保存
' t, o. q* r$ c6 {curve_val(1,:)=x_uni';
  O% ~/ b3 u- a! q5 kcurve_val(2,:)=y_uni;
& K1 N7 a% b' Y%% //对提取出的数据进行拟合(按实际情况进行修改)$ v+ V* H6 E$ X6 G( y2 {7 K
[p,s]=polyfit(curve_val(1,:),curve_val(2,:),4);  %//多项式拟合(为避免龙格库塔,多项式拟合阶数不宜太高)
* }  X+ s- W' H[y_fit,DELTA]=polyval(p,x_uni,s);  %//求拟合后多项式在x_uni对应的y_fit值- ^3 U7 o8 Y5 T7 G: x4 l
figure,plot(x_uni,y_fit),title('拟合后的曲线')) o2 l9 G" l; X. C" |" o9 D
axis([min_x,max_x,min_y,max_y])  %//根据输入设置坐标范围9 ~1 U, R& g! U6 t2 g7 P

该用户从未签到

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

本版积分规则

关闭

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

EDA365公众号

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

GMT+8, 2025-10-31 19:46 , Processed in 0.171875 second(s), 26 queries , Gzip On.

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

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

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