forked from tszalay/PoreView
-
Notifications
You must be signed in to change notification settings - Fork 0
/
find_events.m
103 lines (81 loc) · 3.11 KB
/
find_events.m
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
function DNAevent = find_events(cf)
%FIND_EVENTS Finds events and returns an array of DNAevent structs
% events = find_events(cf)
% This is just an example of how to use CrampFit.
% Alternatively, instead of passing a CrampFit instance to the code,
% you could just initialize one here:
% cf = CrampFit(filename);
% and then add any virtual signals you need, set up the panels, and
% then start finding events on your merry way...
DNAevent = [];
thresh = 0.07;
% use the signal specified in the second panel of CrampFit for the
% event finding. you'd probably just hardcode this in your own code.
sig = cf.psigs(2).sigs;
% loop through entire file, a bit at a time
curind = 0;
while 1
% find next data exceeding threshold, stepping current index
curind = cf.data.findNext(@(d) d(:,sig) > thresh, curind);
% if we didn't find any, we're done with the file
if curind < 0
break
end
imin = curind;
% find the end of the event
imax = cf.data.findNext(@(d) d(:,sig) < 0.75*thresh,curind);
% make sure we have an end for the event
if imax < 0
break
end
% shift event by one sample in each directon to get whole event
imin = imin-1;
imax = imax+1;
% and move curind too
curind = imax;
% event? maybe event?
ts = cf.data.si*[imin imax];
% time range to view, extended past the event a bit
viewt = [ts(1)-0.001, ts(2)+0.001];
% we've decided we have a possibly good event, then
dna = [];
% store the data we want, including times
% note that we're only grabbing the signal we're analyzing
dna.data = cf.data.get(imin:imax,[1 sig]);
% and the start and end times for the event
dna.tstart = ts(1);
dna.tend = ts(2);
% the average current blockage
dna.blockage = abs(mean(dna.data(:,2)));
% now query on-screen, see what we think
% first, zoom in
cf.setView(viewt);
% then, draw some stuff
h = cf.getAxes(2);
plot(h, [viewt(1) ts(1) ts(1) ts(2) ts(2) viewt(2)],...
[0 0 dna.blockage dna.blockage 0 0],'r');
% this line ignores the stuff you drew, in case you're wondering
cf.autoscaleY();
% do this just for kicks
cf.setCursors(ts);
% if you want it to run automatically, this would be the part you
% would change, or something
k = cf.waitKey();
% clear, and force a redraw. this way, you don't accidentally press
% keys twice because you don't know if it's thinking or not.
cf.clearAxes();
pause(0.01);
% handle key input
if (k == 'q')
return
elseif (k ~= 'y')
continue
end
% stupid matlab structs bug
if isempty(DNAevent)
DNAevent = dna;
else
DNAevent(end+1) = dna;
end
end
end