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

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

[复制链接]

该用户从未签到

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

EDA365欢迎您登录!

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

x
1 A* J- ~6 \! h5 C8 m+ x
目录5 c1 t2 c' y& a, t3 K; b
0.引言
$ D9 G# T( u- k1.思路详解与分析: r# A4 t8 l9 B" F! b/ P
2.MATLAB程序' S# `2 o  v3 K
5 H+ S) J2 |- F' s7 ~
0.引言" S0 L/ T7 y( S
  在读文献的时,经常遇到这样的情况:文章里提出的方法好有趣啊,好想拿文中用的数据来试试看看能不能得到相近的结果,可是文中只有根据原始数据绘制的曲线图,没有数据。如下图所示。
6 s+ a1 {, M7 ~7 s: q* l) R0 O5 l$ q8 |2 E' w# P. J2 O( K0 m9 H

4 g) x. G' [7 N8 t. G4 R6 u7 y, T- P( K
  此时,如果能从文中把这幅图截取下来,输入到一个函数中去,最后能返回从图片中提取到的曲线的坐标数据,岂不美哉。这便是本文的工作。: W1 ~. r, F+ r) e
  k! ]* S" k" R* m8 V1 g. i
1.思路详解与分析
3 \9 S2 l6 ]2 m/ X1.1准备待提取数据的曲线图片) w7 X$ @: v& U2 B/ J) U% i
0 M) E: G% {* m7 l0 I5 ^: M# q
将待提取数据的曲线的图片(如.jpg格式图片),利用 imread 输入到 matlab 中。3 Q# h' u% m) x( e/ Y! O+ w. ]1 ]

* Z$ a7 G; y( @+ z1.2曲线图片预处理与数据转换
3 ?; G, d+ n2 q/ W* u& V* k4 T& u4 \0 \/ S" `& ~" B
曲线图片预处理步骤的主要工作包含如下:6 Q- \+ g( N$ V6 U

2 m$ `/ k9 b/ j( j(1)图像二值化; D8 M8 Z; B( |: @& n7 E
8 C: N+ Q* R( Z+ x+ e3 f# u
将输入图像进行二值化处理,但分割得到的结果并不全为数据,其中可能还包括坐标轴等干扰点需要去除。
7 J8 A# d& N) S# G% R9 D+ t) {: U- J5 z5 r( i/ B( {6 r
(2)获取从图片像素到曲线坐标的定标数据
7 X2 x% T; M, S4 d, S4 o; S8 [  N# v8 [4 Y) m# e- U; |. {
首先,通过ginput()手动从图片中提取到两个像素点,这两个点分别为曲线坐标框的左上角和右下角。# F' y% B( ~% F9 T& s! Y

$ B$ p6 }) f  ~( j7 J/ q 9 U3 B! u( C6 `, M% U

% _5 S6 D: C8 ?2 V$ v: f此时,便获得了曲线在图片上的像素范围
' g; ]- @8 P  ?- b
. h0 ?% ]9 p' L0 u9 i; R  t: J[x_index_min, x_index_max] & [y_index_min, y_index_max]  R3 C3 _9 T$ u& g$ d, o% n" s% [( q" `
4 |8 `3 Q4 I+ V! T" u; q( [
然后,手动输入实际曲线的数据坐标范围 [x_min, x_max] & [ymin, y_max_]. 如下所示。
& x& ~. h3 z' H  [
6 D: l9 }+ Y9 i$ [6 Q4 Y- o+ i : q  V4 ~) d! r; z& j

4 c, P, \/ u7 s( y1 [此时,一方面得到了像素坐标,一方面得到了实际坐标。接下来便利用这对数据,将图片中全部的像素坐标转换到实际坐标。
1 Y: k6 z3 ~  {( w  v& @: g3 u) L! O
最终,得到了由图片提取到的数据散点图,如下:  d3 A- l% d( ~! I9 C" H+ {1 S

5 k& L6 m  J0 R2 K" @
$ b) B7 n" p$ r" R, g! L9 A  c可以看到,此时得到的结果,虽然曲线与所需要的相近,但曲线外的部分,如坐标轴框、坐标轴刻度等也被转换成了数据,还需要进一步的处理。# f" A! R. v" H- f: V9 [

1 k- `9 ]+ P% g4 [, ~7 B1.3数据的进一步处理并得到最终曲线
9 M$ Y/ D' U& z4 b9 a, R. k- A9 L9 j! a. Y
这一步的主要工作是在数据中除去曲线之外的部分(包括坐标轴框、坐标轴刻度等);以及解决一个x数据对应多个y数据的情况。
9 X( L$ ~# o; @; ?) o# g
) K  j' I- T" ~  N/ K) q2 f" P) v/ ]& g显然,坐标轴在整幅数据中,均处于边界位置,因此,很容易想到的一种方法是,设定阈值,将距离边界较近的数据直接删除。这里,设定了两个阈值,一个用来限定x方向上的数据,一个用来限定y方向上的数据。比如设定:rate_x=0.08; rate_y=0.05; 意思是阈值设定为曲线最前端8%和最后段8%的数据,曲线最顶端5%和最底端5%的数据。, t  ~9 Z# ^0 Z# L' Z( g: J

1 o4 s, n# C# T, g: G( a' A, }进一步的,对于提取到的数据图,大多数情况一个x会对应若干y,因为数据是由图片转化而来,而图片的分辨率有限,一个实际数据点会用多个像素来表示。解决此问题的中心思想是将同一个x对应的若干个y取均值,但不能直接求均值,因为还有很多y是噪声(如坐标轴线、由图片噪声带入的干扰点等)需要先去除,在第一个问题中,通过限定y的范围,已经在一定程度上除去了部分干扰,在此基础上,我们求取一个x对应的这组y值的均值mean与标准差std,当某些y值位于[mean-std , mean+std]之外,则认为这些y值波动太大,将它们删去。' |. f7 a1 q* z2 ~, n0 Q7 a
1 L+ F+ i: }0 m  N7 q6 U$ J( k
到这里,我们就将数据的处理部分基本完成了,我们将处理后的数据再次绘制成曲线,便可以得到如下/ n$ z! O) s; o0 M( m* S  M0 v

4 b1 L5 {% H) z! T- ^6 t. X
! h2 q6 \: A: M8 `. g  _4 p# ~- a
对比处理之前的数据,由于限定了范围,因此曲线图片中带来的坐标框等内容转化而来的数据已经被删去。
* S" m0 v1 b/ ]5 f6 e
4 N! l& I7 r4 H9 P: e- h将需要提取坐标的曲线图片,和我们提取并处理后的数据绘制的曲线,放在一起对比如下:
9 G% _! B( t3 p: k0 b+ k3 A' I' o& G; _; I' X3 d; z/ a  X

% ~* \! v+ }) \
7 t% p3 a# |( G/ ^  t2 F可以看到,与原曲线图片相比,提取到的数据曲线相似度能达到较高要求。但进一步观察会发现,右图曲线较左图而言,高频分量有一定的减少(即右图曲线更平滑),原因在于数据处理时,对同一个x对应的这组y值进行了均值处理,则在图像上近似反应为均值滤波,从而使得提取到的数据绘制成的曲线的高频分量被滤除。
4 T8 W0 P9 P8 f( u$ Y  H7 y5 z4 l$ I$ L. U
最后,将提取到的最终数据,保存起来如下:
7 B% Q. A/ A8 E; ~; I8 [2 h" {. k# ]
, o! S  G. F+ d' B
( w& L4 G+ L3 m
1.4进一步的讨论——曲线拟合
6 }& h- g/ U% \. ]7 W7 f) e6 ?# }9 V8 |. _' T7 k  d
通过对图片中曲线的数据提取,可以得到数值上的答案,这会带来进一步的思考,即能否得到这些数据的解析表达式。很容易想到,利用最小二乘法来拟合这些数据,这便涉及到了曲线的拟合。(插值与拟合可以这么理解:对于数据点集,若均落在曲线上,则该曲线为插值曲线,否则为拟合曲线)
' ?  n8 X; `: W1 I  y0 I& ^& ?2 i# F1 O6 r& m( p
对于一些简单的曲线图片(如下),可以考虑用泰勒级数来近似,即多项式拟合。2 N  O0 x7 H  H
2 h/ M# }! `8 y( B# b8 @6 d
1 w0 Y7 W# o0 Y# y  A; V- p) N* Q
, c1 c) g, X5 F  [( k  y" y
数据提取并拟合的结果展示如下
1 a1 H$ B- M+ n" `1 A
3 x1 z$ r; a5 |  X4 ] , A7 ]6 k: \6 {& m, C

: t; y& \8 @6 s% B8 o同时还能得到拟合多项式的系数! T% y4 O) }8 d3 q/ w

$ M4 J/ G6 w# M" V
& T4 s+ D8 p( w; {
! v; \( i- v% f1 [4 c. j' a; P从而得到该曲线的多项式(这里采用四阶多项式)表达式为:
* `9 m. L+ T1 w# w+ k4 r
- p$ T/ d: q4 s9 C' F0 a
  h: c4 f( K0 }- k1 d/ C: ?- t, f! ^/ k0 y  b) X
理论上,泰勒级数可以分解任何函数,但实际上,多项式拟合的次数太高,会出现龙格库塔现象,即摆尾现象。因此,多项式拟合的阶数不易过高,一般低于5阶。对于本文最开始的那幅曲线而言,仅5阶的泰勒级数就显得力不从心了,因此,对于这种存在波动剧烈的函数,可以考虑用傅里叶级数进行拟合,或者如果能提供先验知识,可以直接用先验表达式进行拟合。
3 Q8 o2 l# p, z, ]# K, n2 p! G& e" |9 ~% d# W
在MATLAB中,提供了cftool工具箱,其提供了拟合与插值的GUI,使用非常方便,直接在命令窗输入cftool即可调用,cftool界面如下所示,其具体使用方法不在此赘述。4 s' C6 a5 @; h6 j: N

& ~( v2 U( T- T9 ?( c4 L
8 g1 q& m% W2 ]/ W( J$ Y7 ]$ S- X' @5 A7 N' v* [) I0 q+ T
2.MATLAB程序
+ s; X5 C) C2 T- r; T: }" V2 YMATLAB源代码如下所示,和以往的风格一样,提供了详细的注释
3 y, x; C/ V; l# {* [* [3 k4 y* k: ]+ l. D; d6 }
% //提取图片中的曲线数据$ k- b" t4 a# ~
clear,clc,close all
/ @4 H9 k2 h& Y1 l%% //图片与曲线间的定标5 E/ j1 Y* r( }. V* `2 Q$ V) G
im=imread('tu1.jpg'); %//读入图片(替换成需要提取曲线的图片)4 j1 ^: O8 g7 o: s$ L7 f
im=rgb2gray(im); %//灰度变化
7 n  q+ @( t$ Hthresh = graythresh(im); %//二值化阈值
; l+ c; ?$ d8 uim=im2bw(im,thresh); %//二值化5 M$ z" G4 O9 E2 F, D
set(0,'defaultfigurecolor','w')" G0 h% b" G" L0 g) f+ }- U' O
imshow(im) %//显示图片$ c# r4 g; s: t: O/ p% Z  k
[y,x]=find(im==0); %//找出图形中的“黑点”的坐标。该坐标是一维数据。
/ B, c( h! k4 H. ^y=max(y)-y; %//将屏幕坐标转换为右手系笛卡尔坐标
% A2 S* |" D. o) G/ a# i4 w7 ^! J  xy=fliplr(y); %//fliplr()——左右翻转数组
; R; a6 z4 s0 Gplot(x,y,'r.','Markersize', 2);8 z+ f# q5 f5 H3 M! h
disp('请在Figrure中先后点击实际坐标框的两个顶点(左上点和右下点),即A、B两点. ');
) Q0 }; W$ c1 X3 }2 l[Xx,Yy]=ginput(2); %//Xx,Yy——指实际坐标框的两个顶点. w0 m3 x. _0 g1 ~
min_x=input('最小的x值');  %//输入x轴最小值8 q2 M. O' D% `
max_x=input('最大的x值');  %//输入x轴最大值
0 t( J: E6 a: b: o: L+ fmin_y=input('最小的y值');  %//输入y轴最小值7 W# {& I: w# ^
max_y=input('最大的y值');  %//输入y轴最大值
( S! _( Y& S/ n; n- O6 zx=(x-Xx(1))*(max_x-min_x)/(Xx(2)-Xx(1))+min_x;
8 G9 o8 X5 o7 c1 cy=(y-Yy(1))*(min_y-max_y)/(Yy(2)-Yy(1))+max_y;
  O' l/ d& v. M6 v5 n. B, Nplot(x,y,'r.','Markersize', 2);
4 R: h: U; G" H8 ]% Faxis([min_x,max_x,min_y,max_y])  %//根据输入设置坐标范围
8 H% A+ ]6 ~' D. E5 Rtitle('由原图片得到的未处理散点图')
- i- O" k* c6 q/ L9 u%% //将散点转换为可用的曲线
/ z" l# `3 e0 k) {7 `. u%//需处理的问题与解决思路
$ f) i. Y- }, j. t' y( r%//(1)散点图中可能一个x对应好几个y <---> 保留mean()-std()到mean()+std()之间的y值 并取平均处理
6 ^" s/ ~/ M8 B; D6 H0 b& C, a, r%//(2)曲线的最前端和最后段干扰较大 <---> 去掉曲线整体的前(如5%)和后5%8 l* J/ n4 d1 m" i7 V7 I( Z
%//(3)曲线的最顶端和最底段干扰较大 <---> 去掉曲线整体的上10%和下10%
$ \, f+ x  Q6 R) @
: E2 e6 Z" C& K2 L, ^%//参数预设! t! l0 t" Y% _0 d: J7 M
rate_x=0.08;  %//曲线的最前端和最后段删除比例
+ F: i" [( s+ D1 Jrate_y=0.05;  %//曲线的最顶端和最底段删除比例* X/ S( L0 G! G- D8 m

+ ^4 v( B7 C+ K( M7 w7 h2 @[x_uni,index_x_uni]=unique(x);  %//找出有多少个不同的x坐标8 g0 m! f2 t* B. f4 v% |
, M3 T4 ?) h7 E* T( X5 e$ D3 B5 q
x_uni(1:floor(length(x_uni)*rate_x))=[];  %//除去前rate_x(如5%)的x坐标
3 Y, \& a9 D9 v$ rx_uni(floor(length(x_uni)*(1-rate_x)):end)=[];  %//除去后rate_x的x坐标5 Z" l* h7 T8 g2 M6 c
index_x_uni(1:floor(length(index_x_uni)*rate_x))=[];  %//除去前rate_x的x坐标
" c, n( e! @( I' V- Hindex_x_uni(floor(length(index_x_uni)*(1-rate_x)):end)=[];  %//除去后rate_x的x坐标. @1 L" o% ?5 y2 W5 J) B/ @& i
3 w7 {% A8 e/ B+ ]' K# p
[mxu,~]=size(x_uni);! N* x& M& \8 }: R, ~7 N5 U
[mx,~]=size(x);& ^$ f0 ]& C+ V9 j  L  p
for ii=1:mxu& x  n* E4 \+ a6 I4 K' R; a
    if ii==mxu
8 ^4 x9 C9 d' x+ F        ytemp=y(index_x_uni(ii):mx);  G. K! r. S+ L, _7 K% {  C4 f- Z
    else
4 j6 ~: B1 h1 q# N& U( E        ytemp=y(index_x_uni(ii):index_x_uni(ii+1));
0 M! c& m0 R; q# o4 m; d  }3 x    end
) G1 s" W0 A6 H8 N$ ]: n) ~    % //删除方差过大的异常点+ t1 D% x& ^  B" r6 C
    threshold1=mean(ytemp)-std(ytemp);  a2 |# Q6 ?* w; Y' b% S& B
    threshold2=mean(ytemp)+std(ytemp);
/ o) K' M8 i/ J( F    ytemp(find(ytemp<threshold1))=[];  %//删除同一个x对应的一段y中的异常点& W" s/ [' I9 z. D! o
    ytemp(find(ytemp>threshold2))=[];
( V" `) D. t3 l6 c7 A+ o3 a& K2 Z    % //删除距顶端和底端较近的点
. W( @* @! t+ Y& ~    thresholdy=(max_y-min_y)*rate_y;  %//y坐标向阈值# O5 g: E( U2 g
    ytemp(find(ytemp>max_y-thresholdy))=[];  %//删除y轴向距离顶端与底端距离小于rate_y的坐标, K  W: R$ @' I& P
    ytemp(find(ytemp<min_y+thresholdy))=[];  ]4 r4 [0 B& ?0 r! x
    % //剩下的y求均值" Z9 W, L) t- K
    y_uni(ii)=mean(ytemp);+ E* j( u( m& B% Z9 |3 P2 S5 U
end) U/ p( l& n1 A  v$ ?  b2 [
% //此时很多x_uni点处对应的y_uni为空,即NAN,要进一步删去这些空点4 o/ }2 Y/ S( O2 ~& f5 x/ D
x_uni(find(isnan(y_uni)))=[];
* ^( o/ N" H  U4 `# ~y_uni(find(isnan(y_uni)))=[];/ v: j1 }: k7 }% }3 u7 S
% //画图
' [/ M8 C- j% O; L; o7 n+ `6 Qfigure,plot(x_uni,y_uni),title('经处理后得到的扫描曲线')/ R, F  R9 ^9 d8 W( @1 D! {
axis([min_x,max_x,min_y,max_y])  %//根据输入设置坐标范围
+ Z. k% @& t: I1 t' j' _% //将最终提取到的x与y数据保存& @1 I7 S, c8 o5 x9 i( p4 ?
curve_val(1,:)=x_uni';
" j# H- U' ]* S0 r" |curve_val(2,:)=y_uni;* J8 }- T, r9 q3 t
%% //对提取出的数据进行拟合(按实际情况进行修改)
' w* |  H6 R2 K[p,s]=polyfit(curve_val(1,:),curve_val(2,:),4);  %//多项式拟合(为避免龙格库塔,多项式拟合阶数不宜太高)! m. w" t5 i% q) o/ C! V4 A9 [
[y_fit,DELTA]=polyval(p,x_uni,s);  %//求拟合后多项式在x_uni对应的y_fit值
7 B  X4 T! m3 a, ?+ c, Sfigure,plot(x_uni,y_fit),title('拟合后的曲线')8 `7 u: Q5 h/ h* _' z
axis([min_x,max_x,min_y,max_y])  %//根据输入设置坐标范围( }/ B* C& {: O7 A; j/ B- E

该用户从未签到

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

本版积分规则

关闭

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

EDA365公众号

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

GMT+8, 2025-11-24 12:39 , Processed in 0.171875 second(s), 26 queries , Gzip On.

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

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

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