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

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

[复制链接]

该用户从未签到

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

EDA365欢迎您登录!

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

x
* I, [  y2 n! s
目录( r7 H) J7 G  j# }7 L/ k
0.引言
' f9 J. d1 y  r6 s6 Y2 r1.思路详解与分析
! P6 C9 }6 X9 j  z+ x) }7 b. Z8 a; a8 {2.MATLAB程序
. {8 N' C/ v0 N4 { 4 w! G( q8 L! a
0.引言
( K+ b7 ]1 U, I  在读文献的时,经常遇到这样的情况:文章里提出的方法好有趣啊,好想拿文中用的数据来试试看看能不能得到相近的结果,可是文中只有根据原始数据绘制的曲线图,没有数据。如下图所示。
8 c, X7 L. I8 C+ q$ A; X; L& d9 N. l5 v0 ~0 o7 r
1 F) T) |0 a( g
7 @1 X( f% D. K5 P
  此时,如果能从文中把这幅图截取下来,输入到一个函数中去,最后能返回从图片中提取到的曲线的坐标数据,岂不美哉。这便是本文的工作。7 ~9 _/ L5 W$ x7 ~
/ m4 a; P7 d; E7 ?
1.思路详解与分析  v! Q: J8 p' X# t+ J2 E* Q2 j
1.1准备待提取数据的曲线图片; z/ Z4 J& K% d7 _; E1 \

. p9 X$ D5 o5 l  s; [# _- L将待提取数据的曲线的图片(如.jpg格式图片),利用 imread 输入到 matlab 中。1 s8 _. j3 d+ U

. ]  L, E* m7 j1.2曲线图片预处理与数据转换
0 z7 i# [; x! D4 C0 b+ a, ]5 k
6 r1 Z9 _& A3 ]8 h曲线图片预处理步骤的主要工作包含如下:$ ~4 r+ S; w' P0 a) K% X
; @# C* j/ ?; a, T
(1)图像二值化% j( {* D7 L9 H2 z5 e$ X$ {' [

% B4 \& Z' s$ j  m+ d! C1 x* B将输入图像进行二值化处理,但分割得到的结果并不全为数据,其中可能还包括坐标轴等干扰点需要去除。% s) l7 f. s, [+ Q6 {# c* f; ]0 m
* E1 ]- c0 y5 l9 v+ R* r$ R
(2)获取从图片像素到曲线坐标的定标数据# j2 D8 F' ^' {% F- [$ J

) u9 d8 k4 T& w+ }0 P9 L1 _. E首先,通过ginput()手动从图片中提取到两个像素点,这两个点分别为曲线坐标框的左上角和右下角。1 O) T2 I& D2 C  ~
" f% g3 C6 d: m! b$ W0 o* b1 D0 X9 I9 s

& E- T2 K2 G0 V: x5 H' W$ I6 ^1 ^' g3 `4 `* I! l6 N
此时,便获得了曲线在图片上的像素范围
, B" g8 N/ e2 Y# l5 T" _8 q
/ a) f' G3 _1 K, S2 W3 d4 y5 H[x_index_min, x_index_max] & [y_index_min, y_index_max]
' o7 _) w( N: O5 B% m
6 r8 J! F; L, C9 Z然后,手动输入实际曲线的数据坐标范围 [x_min, x_max] & [ymin, y_max_]. 如下所示。
7 B: V5 a% m0 ?9 [6 u  u7 D4 @! h. u' {0 g7 n

  c% s% C' T3 M! o  r  v  i7 S+ D& x9 [$ {
此时,一方面得到了像素坐标,一方面得到了实际坐标。接下来便利用这对数据,将图片中全部的像素坐标转换到实际坐标。
; Q  |" `! d" z* b# B  u" X& X& [+ ?3 y
最终,得到了由图片提取到的数据散点图,如下:
% R) g* {3 b' k  Z6 j, l0 l1 ~) O7 t2 j' ^

" |$ z. d/ {4 U可以看到,此时得到的结果,虽然曲线与所需要的相近,但曲线外的部分,如坐标轴框、坐标轴刻度等也被转换成了数据,还需要进一步的处理。
0 d6 [5 R' l& d' Q
- V& z6 b. h; W. D1.3数据的进一步处理并得到最终曲线2 b5 ~# J8 U  M7 E

9 r! c8 Q' b. F! w这一步的主要工作是在数据中除去曲线之外的部分(包括坐标轴框、坐标轴刻度等);以及解决一个x数据对应多个y数据的情况。
9 `% |0 s8 G: x  v. g
5 h6 h/ J$ `  X! t显然,坐标轴在整幅数据中,均处于边界位置,因此,很容易想到的一种方法是,设定阈值,将距离边界较近的数据直接删除。这里,设定了两个阈值,一个用来限定x方向上的数据,一个用来限定y方向上的数据。比如设定:rate_x=0.08; rate_y=0.05; 意思是阈值设定为曲线最前端8%和最后段8%的数据,曲线最顶端5%和最底端5%的数据。
& Y+ C' b+ w  k/ d9 A. ?: |+ F5 m8 P' Q7 F4 |* O0 \
进一步的,对于提取到的数据图,大多数情况一个x会对应若干y,因为数据是由图片转化而来,而图片的分辨率有限,一个实际数据点会用多个像素来表示。解决此问题的中心思想是将同一个x对应的若干个y取均值,但不能直接求均值,因为还有很多y是噪声(如坐标轴线、由图片噪声带入的干扰点等)需要先去除,在第一个问题中,通过限定y的范围,已经在一定程度上除去了部分干扰,在此基础上,我们求取一个x对应的这组y值的均值mean与标准差std,当某些y值位于[mean-std , mean+std]之外,则认为这些y值波动太大,将它们删去。$ r+ I% ]* M/ [/ ~
/ H3 F3 G& A- X6 R4 _  }
到这里,我们就将数据的处理部分基本完成了,我们将处理后的数据再次绘制成曲线,便可以得到如下
& u, n4 D+ J: O. L! y- V) m' l' J; \( p7 E* `  y
3 y* z  t% p) A5 M( y/ E/ @5 B' F" }
$ f1 Q$ X9 M- w, v( I" b
对比处理之前的数据,由于限定了范围,因此曲线图片中带来的坐标框等内容转化而来的数据已经被删去。) b, m- r. z7 R& G

0 O8 Z3 Z6 J  v将需要提取坐标的曲线图片,和我们提取并处理后的数据绘制的曲线,放在一起对比如下:
2 W5 m" w0 [5 k* M* e
# W. c4 K1 q0 T: x+ S2 Z
* }3 L8 T7 h! w$ ~& f1 h! @5 x) D+ A# q8 Q% ~9 t' V5 |
可以看到,与原曲线图片相比,提取到的数据曲线相似度能达到较高要求。但进一步观察会发现,右图曲线较左图而言,高频分量有一定的减少(即右图曲线更平滑),原因在于数据处理时,对同一个x对应的这组y值进行了均值处理,则在图像上近似反应为均值滤波,从而使得提取到的数据绘制成的曲线的高频分量被滤除。; @9 E$ h9 ?" x. o

4 ?6 r' }4 w1 B1 W  }最后,将提取到的最终数据,保存起来如下:
; N8 G0 q0 o( S2 W) \! m5 E5 F
+ m% X" @) C! E7 b# w+ t + r) ?8 O/ o! M* O. m8 @. D

" o4 Z# }/ ?. m- t1.4进一步的讨论——曲线拟合! \, j* k5 G+ D2 a% Z+ c
! j2 ^1 j& x3 ?+ j) p, h
通过对图片中曲线的数据提取,可以得到数值上的答案,这会带来进一步的思考,即能否得到这些数据的解析表达式。很容易想到,利用最小二乘法来拟合这些数据,这便涉及到了曲线的拟合。(插值与拟合可以这么理解:对于数据点集,若均落在曲线上,则该曲线为插值曲线,否则为拟合曲线)
% b! U" v: t" X! u+ \  O/ h' t- m' i7 w
对于一些简单的曲线图片(如下),可以考虑用泰勒级数来近似,即多项式拟合。
" L  ^3 V+ [! A# J& h" G4 V" p& J/ z
. {2 m! g4 s) \2 x  J8 \" T9 } 2 `  p9 ?2 s( I
$ n3 P- I' n+ S9 t: r: _! i
数据提取并拟合的结果展示如下) l$ [* z* r) P' i- w$ C1 I. d
8 U. g; Y( _* C5 l( H8 U

# |; n$ u7 ]9 c1 ~1 f: a# q0 E; {* i6 R9 v) z4 P
同时还能得到拟合多项式的系数
5 g$ |1 M4 V6 v# V- H9 G5 m% t" v( P: @  D, P! l/ T2 m( z* u5 i
" O. m8 t1 l; `3 Y- C

" w7 e% K1 G* N7 [- e" e4 y从而得到该曲线的多项式(这里采用四阶多项式)表达式为:% G5 r4 U* t) \
0 T/ T, T8 l' t& n% Y7 @
) r! z- a6 Z. X( _% N2 M4 r3 \

' o- s$ w; ^5 |, h4 t7 z6 T理论上,泰勒级数可以分解任何函数,但实际上,多项式拟合的次数太高,会出现龙格库塔现象,即摆尾现象。因此,多项式拟合的阶数不易过高,一般低于5阶。对于本文最开始的那幅曲线而言,仅5阶的泰勒级数就显得力不从心了,因此,对于这种存在波动剧烈的函数,可以考虑用傅里叶级数进行拟合,或者如果能提供先验知识,可以直接用先验表达式进行拟合。
, ?) A" F" P/ `( }. J6 S, q
; C! P, {( Y7 F在MATLAB中,提供了cftool工具箱,其提供了拟合与插值的GUI,使用非常方便,直接在命令窗输入cftool即可调用,cftool界面如下所示,其具体使用方法不在此赘述。/ l; v/ m# e! q: O. F1 u
; V; X1 o2 r6 }( `% n
0 ?5 \5 ]7 w+ S4 T/ N

/ Q7 X) X# s& c( _2.MATLAB程序" J! L$ R8 R4 s% S, _; d( N7 `( @/ I
MATLAB源代码如下所示,和以往的风格一样,提供了详细的注释
7 a! Z- ~: {# ~- ~, f; i! b; G( E/ e' U
% //提取图片中的曲线数据
4 a  [1 p& C1 q; ?2 Y7 tclear,clc,close all
' d4 E/ `8 u8 R3 N%% //图片与曲线间的定标
% q- v- g- R( jim=imread('tu1.jpg'); %//读入图片(替换成需要提取曲线的图片)
6 E1 J' m  N  O/ R& dim=rgb2gray(im); %//灰度变化
. @* H" O; |1 s" J- L4 V: Kthresh = graythresh(im); %//二值化阈值
% Z) G: D+ b4 Zim=im2bw(im,thresh); %//二值化( W, h: Z8 s/ U
set(0,'defaultfigurecolor','w')5 [3 X& z' {/ y. B9 |, b
imshow(im) %//显示图片
3 N- q) F( `, ^/ _9 C9 D  c[y,x]=find(im==0); %//找出图形中的“黑点”的坐标。该坐标是一维数据。
" N3 ~  R2 b# o9 b0 W# {y=max(y)-y; %//将屏幕坐标转换为右手系笛卡尔坐标
7 d; Q: I5 V3 T+ o! [+ K3 }( S4 qy=fliplr(y); %//fliplr()——左右翻转数组+ s+ n) v4 p" D7 C1 L
plot(x,y,'r.','Markersize', 2);
1 X9 @, R& B" {$ R2 w; T2 d- hdisp('请在Figrure中先后点击实际坐标框的两个顶点(左上点和右下点),即A、B两点. ');
- r5 q/ z( }$ F[Xx,Yy]=ginput(2); %//Xx,Yy——指实际坐标框的两个顶点- k, O: h3 a  V0 |& N
min_x=input('最小的x值');  %//输入x轴最小值
, Z2 P$ W" I5 W* M/ kmax_x=input('最大的x值');  %//输入x轴最大值
( N' M( f9 }0 ^0 wmin_y=input('最小的y值');  %//输入y轴最小值
  [2 T: P( |2 S, B- w8 ~1 omax_y=input('最大的y值');  %//输入y轴最大值5 t/ X& W7 |/ t& Y7 U, `% c" U
x=(x-Xx(1))*(max_x-min_x)/(Xx(2)-Xx(1))+min_x;
2 H4 W, o3 Z) \& m- n7 by=(y-Yy(1))*(min_y-max_y)/(Yy(2)-Yy(1))+max_y;4 f3 Y0 L+ F  }2 q
plot(x,y,'r.','Markersize', 2);
' v# x0 M8 p# P1 q: baxis([min_x,max_x,min_y,max_y])  %//根据输入设置坐标范围6 u1 ?! U- [# N+ ^- g3 K
title('由原图片得到的未处理散点图')% W5 o6 D: ^! K$ r. H: [
%% //将散点转换为可用的曲线% V3 L9 C; [  z
%//需处理的问题与解决思路/ L9 C9 V8 {& ~$ m  G
%//(1)散点图中可能一个x对应好几个y <---> 保留mean()-std()到mean()+std()之间的y值 并取平均处理0 F7 A6 r, N. T/ k" A
%//(2)曲线的最前端和最后段干扰较大 <---> 去掉曲线整体的前(如5%)和后5%4 V0 f# a# |( K
%//(3)曲线的最顶端和最底段干扰较大 <---> 去掉曲线整体的上10%和下10%
1 c9 \+ P9 k+ _) x' c% _2 L" d
" n) ?  U$ _. l, X7 u%//参数预设! k! ^: {# v* }7 P  L
rate_x=0.08;  %//曲线的最前端和最后段删除比例% m+ [% B# Q/ n8 ]& x: F
rate_y=0.05;  %//曲线的最顶端和最底段删除比例
+ D" }* L2 }, U( t& Y7 I6 F& S+ W
[x_uni,index_x_uni]=unique(x);  %//找出有多少个不同的x坐标; `3 d$ J6 H! ~8 ^

% }: v) g6 X; Nx_uni(1:floor(length(x_uni)*rate_x))=[];  %//除去前rate_x(如5%)的x坐标
# [5 `" Z4 N9 px_uni(floor(length(x_uni)*(1-rate_x)):end)=[];  %//除去后rate_x的x坐标0 e( N+ w  n- e# O
index_x_uni(1:floor(length(index_x_uni)*rate_x))=[];  %//除去前rate_x的x坐标
' I  j) t2 s" _- y3 x; d/ Nindex_x_uni(floor(length(index_x_uni)*(1-rate_x)):end)=[];  %//除去后rate_x的x坐标- b0 ]1 [- N+ r/ c# u

' R$ i2 v6 t7 k. S  `[mxu,~]=size(x_uni);
5 v7 N- N4 t" L4 ]$ o[mx,~]=size(x);- d- B; T/ Y. |7 g$ g
for ii=1:mxu6 j5 L3 f9 y; V* x  T% [) E
    if ii==mxu
$ D: ]5 |' u3 Z) M* z$ s; w        ytemp=y(index_x_uni(ii):mx);
4 h; T8 L' \/ @8 l+ A& v0 ^% J% I. ]    else
  Q8 M$ w1 d9 p        ytemp=y(index_x_uni(ii):index_x_uni(ii+1));, u; z" [( u6 M! j7 \9 Q- Y% i
    end
( ?' g; b4 n6 }# t% R* z0 ?    % //删除方差过大的异常点6 f7 }# t2 |6 T
    threshold1=mean(ytemp)-std(ytemp);
% ]3 o) _6 s3 d: c2 M' D* u$ V/ ^    threshold2=mean(ytemp)+std(ytemp);) G* i, @/ v; v  S0 ?4 B
    ytemp(find(ytemp<threshold1))=[];  %//删除同一个x对应的一段y中的异常点
  c) s! N8 O' P! F8 H! |    ytemp(find(ytemp>threshold2))=[];, b  h; R7 C1 [+ c3 h
    % //删除距顶端和底端较近的点
+ d( l# c1 Z& v- z# ?    thresholdy=(max_y-min_y)*rate_y;  %//y坐标向阈值
5 Q( J  `8 P- G( g1 A; O* u    ytemp(find(ytemp>max_y-thresholdy))=[];  %//删除y轴向距离顶端与底端距离小于rate_y的坐标3 L3 a' y2 Y+ ?3 h4 y' M
    ytemp(find(ytemp<min_y+thresholdy))=[];
* @* p6 R0 w% ~# h- T    % //剩下的y求均值0 {: K. G" v0 U
    y_uni(ii)=mean(ytemp);
& L! l1 h% b6 P6 q) bend
& \  I( M; e- A' R' h) D% //此时很多x_uni点处对应的y_uni为空,即NAN,要进一步删去这些空点" O5 p6 `) a8 ~% A: L: m5 r
x_uni(find(isnan(y_uni)))=[];
% S) h1 [+ k  u: {$ vy_uni(find(isnan(y_uni)))=[];# K# }* Z/ ]6 |2 H+ q& T
% //画图
: P9 O) W( }0 h3 l; \" U6 f/ kfigure,plot(x_uni,y_uni),title('经处理后得到的扫描曲线')) C8 M' ]4 E2 j- `: ?2 U
axis([min_x,max_x,min_y,max_y])  %//根据输入设置坐标范围6 D/ O& z; n% N/ t, B; A8 D. l1 l- X
% //将最终提取到的x与y数据保存0 b6 p3 V  O$ F/ G% S- g& p
curve_val(1,:)=x_uni';
! s/ S8 ~- r2 t, D% fcurve_val(2,:)=y_uni;9 e2 L5 n2 K3 m* W
%% //对提取出的数据进行拟合(按实际情况进行修改)
$ f4 m  |- q! }7 I5 L" x5 z% _[p,s]=polyfit(curve_val(1,:),curve_val(2,:),4);  %//多项式拟合(为避免龙格库塔,多项式拟合阶数不宜太高). }$ ~0 K) ]! q0 _. J
[y_fit,DELTA]=polyval(p,x_uni,s);  %//求拟合后多项式在x_uni对应的y_fit值
7 g0 M6 G/ k6 m4 ~3 sfigure,plot(x_uni,y_fit),title('拟合后的曲线')
$ m5 q0 t) G- ~! V$ P) r/ E/ Laxis([min_x,max_x,min_y,max_y])  %//根据输入设置坐标范围, S4 x7 G5 ^( |9 O' N/ ]: X

该用户从未签到

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

本版积分规则

关闭

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

EDA365公众号

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

GMT+8, 2025-8-5 03:04 , Processed in 0.125000 second(s), 26 queries , Gzip On.

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

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

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