-
Notifications
You must be signed in to change notification settings - Fork 1
/
PreprocessSpectrum.m
79 lines (68 loc) · 2.21 KB
/
PreprocessSpectrum.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
function [w_out,r_out] = PreprocessSpectrum(w,r)
% Preprocess a given spectrum to be used in the clustering algorithm.
% Returns the preprocessed spectrum, or an empty if preprocessing
% conditions do not hold.
%
% Arguments:
% 'w':
% A column vector with the wavelengths, in units um (micrometer).
% 'r':
% A column vector of the same size as 'w', with reflectances.
%
% Returns:
% 'w_out':
% A column vector with the wavelengths 0.5um to 2.43um, in 0.02um
% steps.
% 'r_out':
% A column vector the same size as 'w_out' with the reflecatances
% after preprocessing. It is obtained by smoothing 'r' and
% interpolating it linearly to the wavelengths in 'w_out'.
%
% If the range of wavelengths does not include the range 0.5um to 2.43um,
% or there is a gap larger than 0.1um in that range, the spectrum cannot be
% processed and the output is empty vectors.
start_wavelength = 0.5; % um
stop_wavelength = 2.43; % um
step_wavelength = 0.02; % um
max_wavelength_gap = 0.1; % um
[w,inds] = sort(w);
r = r(inds);
% remove reflectances that clearly indicate an error in measurement
w = w(r>0.2);
r = r(r>0.2);
if numel(w) < 3
w_out = [];
r_out = [];
return;
end
% if the interpolation range is not contained in the spectrum, can't
% proceed
if (w(1) > start_wavelength) || (w(end) < stop_wavelength)
w_out = [];
r_out = [];
return;
end
% for wavelengths with multiple values, take the average
[w,~,inds] = unique(w);
r = accumarray(inds,r)./accumarray(inds,1);
% if there is a wavelength gap too big, can't proceed
largest_gap = max(diff(w((w>=start_wavelength) & (w<=stop_wavelength))));
if (~isempty(largest_gap) && (largest_gap>max_wavelength_gap))
w_out = [];
r_out = [];
return;
end
% smooth prior to interpolation, using a quadratic kernel with width of
% the step size of interpolation
r_smoothed = r;
for ii=1:numel(w)
inds = find(abs(w-w(ii))<step_wavelength/2);
weights = (step_wavelength/2-abs(w(inds)-w(ii))).^2;
weights = weights/sum(weights);
r_smoothed(ii) = sum(r(inds).*weights);
end
r = r_smoothed;
% interpolate
w_out = [start_wavelength:step_wavelength:stop_wavelength];
r_out = interp1(w,r,w_out,'linear');
end