-
Notifications
You must be signed in to change notification settings - Fork 1
/
Copy pathKalman_Stack_Filter.m
121 lines (117 loc) · 4.12 KB
/
Kalman_Stack_Filter.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
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
function imageStack=Kalman_Stack_Filter(imageStack,gain,percentvar)
% function imageStack=Kalman_Stack_Filter(imageStack,percentvar,gain)
%
% Purpose
% Implements a predictive Kalman-like filter in the time domain of the image
% stack. Algorithm taken from Java code by C.P. Mauer.
% http://rsb.info.nih.gov/ij/plugins/kalman.html
%
% Inputs
% imageStack - a 3d matrix comprising of a noisy image sequence. Time is
% the 3rd dimension.
% gain - the strength of the filter [0 to 1]. Larger gain values means more
% aggressive filtering in time so a smoother function with a lower
% peak. Gain values above 0.5 will weight the predicted value of the
% pixel higher than the observed value.
% percentvar - the initial estimate for the noise [0 to 1]. Doesn't have
% much of an effect on the algorithm.
%
% Output
% imageStack - the filtered image stack
%
% Note:
% The time series will look noisy at first then become smoother as the
% filter accumulates evidence.
%
% Rob Campbell, August 2009
% Process input arguments
if nargin<2, gain=0.5; end
if nargin<3, percentvar = 0.05; end
if gain>1.0||gain<0.0
gain = 0.8;
end
if percentvar>1.0 || percentvar<0.0
percentvar = 0.05;
end
%Copy the last frame onto the end so that we filter the whole way
%through
imageStack(:,:,end+1)=imageStack(:,:,end);
%Set up variables
width = size(imageStack,1);
height = size(imageStack,2);
stacksize = size(imageStack,3);
tmp=ones(width,height);
%Set up priors
predicted = imageStack(:,:,1);
predictedvar = tmp*percentvar;
noisevar=predictedvar;
%Now conduct the Kalman-like filtering on the image stack
for i=2:stacksize-1
stackslice = imageStack(:,:,i+1);
observed = stackslice;
Kalman = predictedvar ./ (predictedvar+noisevar);
corrected = gain*predicted + (1.0-gain)*observed + Kalman.*(observed-predicted);
correctedvar = predictedvar.*(tmp - Kalman);
predictedvar = correctedvar;
predicted = corrected;
imageStack(:,:,i)=corrected;
end
imageStack(:,:,end)=[];function imageStack=Kalman_Stack_Filter(imageStack,gain,percentvar)
% function imageStack=Kalman_Stack_Filter(imageStack,percentvar,gain)
%
% Purpose
% Implements a predictive Kalman-like filter in the time domain of the image
% stack. Algorithm taken from Java code by C.P. Mauer.
% http://rsb.info.nih.gov/ij/plugins/kalman.html
%
% Inputs
% imageStack - a 3d matrix comprising of a noisy image sequence. Time is
% the 3rd dimension.
% gain - the strength of the filter [0 to 1]. Larger gain values means more
% aggressive filtering in time so a smoother function with a lower
% peak. Gain values above 0.5 will weight the predicted value of the
% pixel higher than the observed value.
% percentvar - the initial estimate for the noise [0 to 1]. Doesn't have
% much of an effect on the algorithm.
%
% Output
% imageStack - the filtered image stack
%
% Note:
% The time series will look noisy at first then become smoother as the
% filter accumulates evidence.
%
% Rob Campbell, August 2009
% Process input arguments
if nargin<2, gain=0.5; end
if nargin<3, percentvar = 0.05; end
if gain>1.0||gain<0.0
gain = 0.8;
end
if percentvar>1.0 || percentvar<0.0
percentvar = 0.05;
end
%Copy the last frame onto the end so that we filter the whole way
%through
imageStack(:,:,end+1)=imageStack(:,:,end);
%Set up variables
width = size(imageStack,1);
height = size(imageStack,2);
stacksize = size(imageStack,3);
tmp=ones(width,height);
%Set up priors
predicted = imageStack(:,:,1);
predictedvar = tmp*percentvar;
noisevar=predictedvar;
%Now conduct the Kalman-like filtering on the image stack
for i=2:stacksize-1
stackslice = imageStack(:,:,i+1);
observed = stackslice;
Kalman = predictedvar ./ (predictedvar+noisevar);
corrected = gain*predicted + (1.0-gain)*observed + Kalman.*(observed-predicted);
correctedvar = predictedvar.*(tmp - Kalman);
predictedvar = correctedvar;
predicted = corrected;
imageStack(:,:,i)=corrected;
end
imageStack(:,:,end)=[];