ImaGIN_AlignSegmentsIter.m
3.8 KB
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
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
function [ alignments, alignedStimulations ] = ImaGIN_AlignSegmentsIter( d, stimulations, segmentHWidth, templateHWidth, corrHWidth, do_plot, iterMax )
if size(d,1) ~= 1
error( 'Invalid d' )
end
fprintf( 1, [ 'alignSegments... ' ] ) ;
alignments = nan( length(stimulations), iterMax ) ;
alignedStimulations = stimulations ;
segments = nan( length(alignedStimulations), 2*segmentHWidth+1 ) ;
for i=1:iterMax
% fprintf( 1, [ 'alignSegments iter ' num2str(i) '...' '\n' ] ) ;
for s=1:length(alignedStimulations)
segments(s,:) = d( alignedStimulations(s) + (-segmentHWidth:segmentHWidth) ) ;
end
[ alignment, ~ ] = alignSegments( segments, templateHWidth, corrHWidth, do_plot ) ;
alignments(:,i) = alignment ;
if isempty(find(alignment,1))
% fprintf( 1, [ 'alignSegments iter ' num2str(i) ' => BREAK' '\n' ] ) ;
break
else
% it happens that aligned segment fall outside full range
% in this case, the segment is not aligned
for s=1:length(alignedStimulations)
if alignedStimulations(s)+alignment(s)+segmentHWidth > length(d)
alignment(s) = 0 ;
end
end
alignedStimulations = alignedStimulations + alignment' ;
end
% if i == iterMax
% error( [ '>> WARNING: alignment not achieved after iterMax = ' num2str(iterMax) ] ) ;
% % fprintf( 1, [ '>> WARNING: alignment not achieved after iterMax = ' num2str(iterMax) '\n' ] ) ;
% end
end
fprintf( 1, [ '(' num2str(i) ' iterations)' '\n' ] ) ;
end
function [ alignment, alignedSegments ] = alignSegments( segments, templateHWidth, corrHWidth, do_plot )
if ~ ( corrHWidth > templateHWidth )
error( 'Invalid corrHWidth' )
end
segmentCount = size(segments,1) ;
segmentWidth = size(segments,2) ;
segmentCenter = round(segmentWidth/2);
templateWidth = 2*templateHWidth+1 ;
corrWidth = 2*corrHWidth+1 ;
corrCount = corrWidth - templateWidth + 1 ;
corr = nan( segmentCount, corrCount ) ;
template = mean( segments(:,segmentCenter+(-templateHWidth:templateHWidth)), 1 ) ;
for s=1:segmentCount
for c=1:corrCount
segStart = segmentCenter - corrHWidth -1 + c ;
segEnd = segStart + templateWidth - 1 ;
r = corrcoef( template, segments(s,segStart:segEnd) ) ;
corr(s,c) = r(2,1) ;
end
end
[ ~, i ] = max( corr,[],2) ;
alignment = segmentCenter - corrHWidth -1 + i + templateHWidth - segmentCenter ;
alignedSegments = applyAlignment( segments, alignment ) ;
if do_plot
plotSegmentsAlignment( segments, alignedSegments, segmentCenter ) ;
end
end
function alignedSegments = applyAlignment( segments, alignment )
alignedSegments = nan( size(segments) ) ;
for s=1:size(segments,1)
if alignment(s) > 0
alignedSegments(s,1:end-alignment(s)) = segments(s,1+alignment(s):end) ;
elseif alignment(s) < 0
alignedSegments(s,1+abs(alignment(s)):end) = segments(s,1:end-abs(alignment(s))) ;
else
alignedSegments(s,:) = segments(s,:);
end
end
end
function plotSegmentsAlignment( segments, alignedSegments, segmentCenter )
figure ;
ax1 = subplot(2,1,1); hold on ;
plot( 1:size(segments,2), segments ) ;
plot( 1:size(segments,2), mean(segments,1), 'color', 'r', 'linewidth', 4 ) ;
line( [ segmentCenter; segmentCenter ], get( gca, 'YLim' )', 'color', 'm', 'linewidth', 2 ) ;
xlim( [ 1 size(segments,2) ] ) ;
title( 'templates' )
ax2 = subplot(2,1,2); hold on ;
plot( 1:size(segments,2), alignedSegments ) ;
plot( 1:size(segments,2), mean(alignedSegments,1), 'color', 'r', 'linewidth', 4 ) ;
line( [ segmentCenter; segmentCenter ], get( gca, 'YLim' )', 'color', 'm', 'linewidth', 2 ) ;
xlim( [ 1 size(segments,2) ] ) ;
title( 'templates aligned' )
linkaxes( [ ax1, ax2 ], 'xy' )
end