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

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

[复制链接]

该用户从未签到

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

EDA365欢迎您登录!

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

x
% x5 i- ], m2 D9 G1 O
目录
6 M0 z' j  k7 E- ?, ?- l0.引言
% L# s8 f) y  `/ |1.思路详解与分析9 E' A) e- G& V9 U+ J6 [) W
2.MATLAB程序0 m! S0 O  b2 }+ _5 q" H
4 |4 I: X3 k) I' ~3 ]# w
0.引言1 L, v7 ^; K; h" V  t/ {0 ?1 x
  在读文献的时,经常遇到这样的情况:文章里提出的方法好有趣啊,好想拿文中用的数据来试试看看能不能得到相近的结果,可是文中只有根据原始数据绘制的曲线图,没有数据。如下图所示。; f* k0 [+ I& `6 ]6 L
$ A: I) ^' S4 H# j/ t; u! Q
0 h& t  N) e' F. m! h  O$ ]$ {

. Q6 _8 k: @5 G0 j: V0 ?  此时,如果能从文中把这幅图截取下来,输入到一个函数中去,最后能返回从图片中提取到的曲线的坐标数据,岂不美哉。这便是本文的工作。
* d$ [+ e5 o) z* x
5 _& d; C9 j+ r! B- e& ^4 Q( U2 o# Q1.思路详解与分析  i' }. z8 n: K* W7 \1 o" [; X# ~
1.1准备待提取数据的曲线图片
' _' v& Q3 o+ s% Q. X6 Z! W% i
- ~" {! B2 i* ?# Z" ]将待提取数据的曲线的图片(如.jpg格式图片),利用 imread 输入到 matlab 中。
' R9 s' |. a9 p( C: D* G( D6 a- t4 r/ E2 O, ]0 Z+ W% @
1.2曲线图片预处理与数据转换
/ y; P6 }/ s1 h: i
, w4 y5 t9 @/ y3 x2 A) E9 t9 c曲线图片预处理步骤的主要工作包含如下:
- l" X: R# M/ N- j; J5 n
8 i" m+ m+ U! W: H(1)图像二值化
$ K- e3 A6 a; H3 U! a+ Q3 T' }* j- k# K7 P5 B: H% e+ Q& j
将输入图像进行二值化处理,但分割得到的结果并不全为数据,其中可能还包括坐标轴等干扰点需要去除。
$ C  C! ~( f" ~6 H# k5 r$ G
/ w% b( h# L& \(2)获取从图片像素到曲线坐标的定标数据
8 y: V1 `/ X; [* l/ P7 a
& V. V1 }* M  [. U+ O& g% f  u首先,通过ginput()手动从图片中提取到两个像素点,这两个点分别为曲线坐标框的左上角和右下角。( M3 l2 e: F2 l% R
) X3 k* s+ e; _! R" B

6 f, ]/ n7 f' D. r
' N' l4 Y$ B; p/ C+ a+ K" w- O此时,便获得了曲线在图片上的像素范围- f5 R0 @' ~& K5 G4 m' O

7 T$ d: S& I5 O; q# b$ y* W/ z( ][x_index_min, x_index_max] & [y_index_min, y_index_max]; l* _/ i2 t2 I/ x7 g

1 e0 u5 U& Q' h% E$ M! E. M1 ?然后,手动输入实际曲线的数据坐标范围 [x_min, x_max] & [ymin, y_max_]. 如下所示。' j# e3 L  }" q3 d3 m* N4 @# z
! Q: @7 x2 [! H; U. O' a
6 C4 |  O  X  e: G/ k3 f

8 I. x- c- D. f% \/ k0 [# R此时,一方面得到了像素坐标,一方面得到了实际坐标。接下来便利用这对数据,将图片中全部的像素坐标转换到实际坐标。+ x% @+ \$ {! Z& O, j* ~

2 W: z3 E/ e8 D5 K0 u最终,得到了由图片提取到的数据散点图,如下:$ U+ P) _" \' m9 a

; Z3 l2 U! z8 H8 ^ 4 [! R8 E% W% Z) M. i
可以看到,此时得到的结果,虽然曲线与所需要的相近,但曲线外的部分,如坐标轴框、坐标轴刻度等也被转换成了数据,还需要进一步的处理。# a2 |1 r# S2 {* A3 t/ z

( [  h. R, [  }$ k' f  v1 Q1.3数据的进一步处理并得到最终曲线
4 p- U, [, ^1 h5 J, }, k( m: N- ~4 A8 q( M
这一步的主要工作是在数据中除去曲线之外的部分(包括坐标轴框、坐标轴刻度等);以及解决一个x数据对应多个y数据的情况。( Z( x* H: ]2 u

( W  F' z( x: E6 x* M显然,坐标轴在整幅数据中,均处于边界位置,因此,很容易想到的一种方法是,设定阈值,将距离边界较近的数据直接删除。这里,设定了两个阈值,一个用来限定x方向上的数据,一个用来限定y方向上的数据。比如设定:rate_x=0.08; rate_y=0.05; 意思是阈值设定为曲线最前端8%和最后段8%的数据,曲线最顶端5%和最底端5%的数据。$ `/ k7 K7 @0 p+ e" t
8 h0 h5 ?0 ~/ k6 N2 y7 W7 A
进一步的,对于提取到的数据图,大多数情况一个x会对应若干y,因为数据是由图片转化而来,而图片的分辨率有限,一个实际数据点会用多个像素来表示。解决此问题的中心思想是将同一个x对应的若干个y取均值,但不能直接求均值,因为还有很多y是噪声(如坐标轴线、由图片噪声带入的干扰点等)需要先去除,在第一个问题中,通过限定y的范围,已经在一定程度上除去了部分干扰,在此基础上,我们求取一个x对应的这组y值的均值mean与标准差std,当某些y值位于[mean-std , mean+std]之外,则认为这些y值波动太大,将它们删去。! {# b1 ?5 I* U7 m6 J+ l' M; D
( W$ x% X- e3 w
到这里,我们就将数据的处理部分基本完成了,我们将处理后的数据再次绘制成曲线,便可以得到如下4 m( X' k! [3 C7 W

' r( G9 J( [  w ; ?5 ^9 i+ X& K/ l

0 B% S) L$ d7 O对比处理之前的数据,由于限定了范围,因此曲线图片中带来的坐标框等内容转化而来的数据已经被删去。' X" @- S0 Y9 t8 z' p- [- M

  z8 Y2 \2 G* R0 T1 b1 }1 |将需要提取坐标的曲线图片,和我们提取并处理后的数据绘制的曲线,放在一起对比如下:
  f6 l: F- _" ^/ @) w" {- G. O9 l% a' `2 q9 H" M: v
: K( r! h$ K' |7 u. L
0 J- ^9 x% f4 d8 f! p! E, Y
可以看到,与原曲线图片相比,提取到的数据曲线相似度能达到较高要求。但进一步观察会发现,右图曲线较左图而言,高频分量有一定的减少(即右图曲线更平滑),原因在于数据处理时,对同一个x对应的这组y值进行了均值处理,则在图像上近似反应为均值滤波,从而使得提取到的数据绘制成的曲线的高频分量被滤除。) @7 Z: [- V" P7 B% y

# `* @9 }5 \; F" H最后,将提取到的最终数据,保存起来如下:
6 O1 B! |5 }  ^
9 |+ n1 H+ G6 x3 s 7 b0 p7 A, @- c! h
) [+ M& E0 \1 _) d6 `
1.4进一步的讨论——曲线拟合* b! u4 u3 R# b  H( Z% H2 y, W
$ k4 Z3 e: v; S7 L1 z6 Y
通过对图片中曲线的数据提取,可以得到数值上的答案,这会带来进一步的思考,即能否得到这些数据的解析表达式。很容易想到,利用最小二乘法来拟合这些数据,这便涉及到了曲线的拟合。(插值与拟合可以这么理解:对于数据点集,若均落在曲线上,则该曲线为插值曲线,否则为拟合曲线)! b3 n, u, {  ^. T" F
  ?0 t: H  R7 b. \; o
对于一些简单的曲线图片(如下),可以考虑用泰勒级数来近似,即多项式拟合。
4 K; Q5 S7 U1 e
$ `) |9 r# U# H# l; z7 A
9 K' M" A" s  K8 A  T; P/ H/ I/ C; a- i2 ~
数据提取并拟合的结果展示如下
% n/ j8 N8 W0 A* T1 c. f9 y. I" Z) m! V6 L/ \0 Q7 L
8 Y3 }3 Z7 f4 u0 F

; e/ T" Z9 s, I, t. e9 K同时还能得到拟合多项式的系数5 q1 c. L; h( n, ?; _9 j9 V
) `' H, e: d8 y7 B6 [; W
' r/ t8 ~7 a( K, x
; }2 I$ v* r0 e9 G) i/ k' J
从而得到该曲线的多项式(这里采用四阶多项式)表达式为:
4 T/ d3 k( u" R" Q; f7 `9 _4 q  ]6 w$ [

* g4 L2 E& r: M' U8 m4 S9 ]9 W% H7 p' t6 [
理论上,泰勒级数可以分解任何函数,但实际上,多项式拟合的次数太高,会出现龙格库塔现象,即摆尾现象。因此,多项式拟合的阶数不易过高,一般低于5阶。对于本文最开始的那幅曲线而言,仅5阶的泰勒级数就显得力不从心了,因此,对于这种存在波动剧烈的函数,可以考虑用傅里叶级数进行拟合,或者如果能提供先验知识,可以直接用先验表达式进行拟合。
" m% h! d5 a9 w' @+ |# B1 ?! C- t5 h" ^) B
在MATLAB中,提供了cftool工具箱,其提供了拟合与插值的GUI,使用非常方便,直接在命令窗输入cftool即可调用,cftool界面如下所示,其具体使用方法不在此赘述。3 U" z1 w2 ]4 N9 Y

! b* p- Y3 y) Q1 a/ Q- w3 ^
0 Y3 R; ~: P- A) d6 G( E- H$ \4 K/ S0 Z+ f$ L) G2 ~2 {
2.MATLAB程序
& M! D& f, \4 z% u* |+ dMATLAB源代码如下所示,和以往的风格一样,提供了详细的注释
# H# |2 @! R1 P4 R. _! u
, R7 R5 d* ]6 B1 ^% //提取图片中的曲线数据7 U- ~" m. M5 q, C* u
clear,clc,close all
: \/ F9 `4 ?' `2 g) G; E%% //图片与曲线间的定标7 n& @; g( v$ m+ g# Z$ ~  t7 H
im=imread('tu1.jpg'); %//读入图片(替换成需要提取曲线的图片)
" S! @9 A% k: n% H  ?im=rgb2gray(im); %//灰度变化- {- b& v5 R6 ]) R% H8 R
thresh = graythresh(im); %//二值化阈值3 K9 \1 p- N# _) I1 }
im=im2bw(im,thresh); %//二值化
  M) d' v8 M' |' k6 s+ t7 iset(0,'defaultfigurecolor','w')) K0 B- M! Z- l& M3 g
imshow(im) %//显示图片$ h, c, R  q' o+ ^" P
[y,x]=find(im==0); %//找出图形中的“黑点”的坐标。该坐标是一维数据。1 ^% w$ M- ~& y6 T, }
y=max(y)-y; %//将屏幕坐标转换为右手系笛卡尔坐标
& g" g. r9 U9 ~2 By=fliplr(y); %//fliplr()——左右翻转数组
2 Q* C) E" \8 K2 R( I$ Y4 ~plot(x,y,'r.','Markersize', 2);4 Q% i( |$ ?# u% F. x! F. \
disp('请在Figrure中先后点击实际坐标框的两个顶点(左上点和右下点),即A、B两点. ');% W( G6 d, [0 V+ _* A: @9 A: a
[Xx,Yy]=ginput(2); %//Xx,Yy——指实际坐标框的两个顶点
; V! D9 m5 p  h2 i& wmin_x=input('最小的x值');  %//输入x轴最小值2 x+ S! W% f- w  K' C
max_x=input('最大的x值');  %//输入x轴最大值  G: I! d! j8 y& i' I
min_y=input('最小的y值');  %//输入y轴最小值1 x. g& O+ X7 H* _8 l
max_y=input('最大的y值');  %//输入y轴最大值+ e- W; Y0 f- d
x=(x-Xx(1))*(max_x-min_x)/(Xx(2)-Xx(1))+min_x;
/ o# x6 e+ a6 H* Ey=(y-Yy(1))*(min_y-max_y)/(Yy(2)-Yy(1))+max_y;5 v( C* |2 K% d3 L  [0 Y, U8 M0 I
plot(x,y,'r.','Markersize', 2);; ]: u) y1 D6 n# A, w; y, {$ ~: P8 U
axis([min_x,max_x,min_y,max_y])  %//根据输入设置坐标范围
5 I& _& S  A% u, stitle('由原图片得到的未处理散点图')
" w. t# [* Z' E7 ~% V; r%% //将散点转换为可用的曲线- w& V+ Q, L2 A3 h
%//需处理的问题与解决思路: y, S2 q2 |# c1 A
%//(1)散点图中可能一个x对应好几个y <---> 保留mean()-std()到mean()+std()之间的y值 并取平均处理3 e  N3 I# s% R+ h. X- k6 n
%//(2)曲线的最前端和最后段干扰较大 <---> 去掉曲线整体的前(如5%)和后5%$ z, @3 t* R5 O" S% {/ X9 [
%//(3)曲线的最顶端和最底段干扰较大 <---> 去掉曲线整体的上10%和下10%
& ^- `" p3 L! n- U+ \, Y% Y" Y5 r" Y; M0 u3 Y6 G
%//参数预设
* x+ M. H1 Z2 }( q; R5 frate_x=0.08;  %//曲线的最前端和最后段删除比例: H( r  H5 U- U7 d0 T( N
rate_y=0.05;  %//曲线的最顶端和最底段删除比例* w: A( C6 m. n) w% h( V! L
0 F% g, z8 m8 P/ S
[x_uni,index_x_uni]=unique(x);  %//找出有多少个不同的x坐标
' I2 v: d8 L  U9 y6 W# @
* ]" ]4 O9 B1 S3 vx_uni(1:floor(length(x_uni)*rate_x))=[];  %//除去前rate_x(如5%)的x坐标
% I: G6 [. z. c5 T- Mx_uni(floor(length(x_uni)*(1-rate_x)):end)=[];  %//除去后rate_x的x坐标1 G3 R9 l8 t( Z
index_x_uni(1:floor(length(index_x_uni)*rate_x))=[];  %//除去前rate_x的x坐标
; |& ^) o9 f: q5 m; C# E* c( Aindex_x_uni(floor(length(index_x_uni)*(1-rate_x)):end)=[];  %//除去后rate_x的x坐标& G- K# m5 u  r! W) a

7 g# H# v: g7 A( o[mxu,~]=size(x_uni);; F7 Q  i1 r+ B( x0 I
[mx,~]=size(x);
4 x4 O$ ~, A3 G/ c6 `6 G! ifor ii=1:mxu
9 ]! G( B- U  H- N    if ii==mxu4 e; |8 w" s9 p% a. p% r
        ytemp=y(index_x_uni(ii):mx);) n& t: R. M5 [' w! I
    else
1 r4 G# W: _2 U% j1 M, P& w; h3 F3 O3 q        ytemp=y(index_x_uni(ii):index_x_uni(ii+1));
6 v! f" \! {9 u8 X) f/ |    end4 Q, C: L) V/ o. R; d# t1 D: k
    % //删除方差过大的异常点
& X4 S2 z& Q  R5 r* [6 k& _    threshold1=mean(ytemp)-std(ytemp);
; a# [8 A/ H) l! r' t1 i0 J* I    threshold2=mean(ytemp)+std(ytemp);# L' a9 x) x  W7 v$ j" Y
    ytemp(find(ytemp<threshold1))=[];  %//删除同一个x对应的一段y中的异常点% i  z" G# r& L
    ytemp(find(ytemp>threshold2))=[];1 |3 s' w: w: y0 D4 L0 O1 B6 _1 Q0 }% S
    % //删除距顶端和底端较近的点
0 F/ p. `( l. c    thresholdy=(max_y-min_y)*rate_y;  %//y坐标向阈值
: g& e/ \- }, q    ytemp(find(ytemp>max_y-thresholdy))=[];  %//删除y轴向距离顶端与底端距离小于rate_y的坐标
) L9 w2 `1 B3 L+ r9 x0 y) T    ytemp(find(ytemp<min_y+thresholdy))=[];
0 i3 j% F/ X5 M1 l$ ~* C5 @6 l    % //剩下的y求均值: a/ c' `. [5 _2 a1 |( B9 N
    y_uni(ii)=mean(ytemp);9 m& q4 x1 w8 q, r! P
end
$ k: r! {- w; `( Y! C9 C% //此时很多x_uni点处对应的y_uni为空,即NAN,要进一步删去这些空点6 O4 N. o- b: l: x" X& f
x_uni(find(isnan(y_uni)))=[];
4 w% O- o6 I, W: O. K+ S6 yy_uni(find(isnan(y_uni)))=[];
2 G/ k8 M1 r2 ]% //画图, z6 B. `% d1 [  x% l! d
figure,plot(x_uni,y_uni),title('经处理后得到的扫描曲线')2 d2 y- `7 r1 c; x4 T8 |
axis([min_x,max_x,min_y,max_y])  %//根据输入设置坐标范围: p. r& z$ n% q6 B
% //将最终提取到的x与y数据保存
7 B5 D  T, Z  H& ~1 zcurve_val(1,:)=x_uni';
; ?/ D/ B$ T. K9 L, \0 Ccurve_val(2,:)=y_uni;% F7 F8 M# U4 ?8 N& \; t
%% //对提取出的数据进行拟合(按实际情况进行修改)
3 J- n  b$ i# p4 U[p,s]=polyfit(curve_val(1,:),curve_val(2,:),4);  %//多项式拟合(为避免龙格库塔,多项式拟合阶数不宜太高)
8 W4 K1 V- Q9 ]6 a: ^' z0 _[y_fit,DELTA]=polyval(p,x_uni,s);  %//求拟合后多项式在x_uni对应的y_fit值
/ `) `2 }. U0 M5 y' lfigure,plot(x_uni,y_fit),title('拟合后的曲线'); }' y. ?3 W, q
axis([min_x,max_x,min_y,max_y])  %//根据输入设置坐标范围
: n% M: T8 y8 i, i2 e

该用户从未签到

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

本版积分规则

关闭

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

EDA365公众号

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

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

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

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

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