%FEATUREDEPTH Plots spectrum in region of an absorption feature and estimates its depth.
% [w, D] = featuredepth(S, w1, w2) finds the local minimum within the
% wavelength region w1 to w2 estimates the depth of minimum relative to
% a straight line (continuum) plotted between w1 and w2.
%
% SEE ALSO
% The definition of feature depth is given in Clark et al.
%
function [featureWavelength, depth] = featuredepth(spectrum, leftWavelength, rightWavelength)
plotWavelengthRange = [leftWavelength - 200, rightWavelength + 200];
plotReflectanceRange = [ 0, 1 ];
plot(spectrum.wavelength, spectrum.data, 'Color', 'b');
axis([plotWavelengthRange plotReflectanceRange]);
xlabel('Wavelength in nm');
ylabel('Reflectance');
if isfield(spectrum, 'name')
title(spectrum.name);
else
title('Specturm of unknown sample');
end
% Mark on the continuum interval points.
leftReflectance = spline(spectrum.wavelength, spectrum.data, leftWavelength);
rightReflectance = spline(spectrum.wavelength, spectrum.data, rightWavelength);
line(leftWavelength, leftReflectance, 'Marker', 's', 'Color', 'r');
line(rightWavelength, rightReflectance, 'Marker', 's', 'Color', 'r');
% Calculate the continuum (which is a straight line).
gradient = (rightReflectance - leftReflectance) / (rightWavelength - leftWavelength);
intercept = (leftReflectance*rightWavelength - rightReflectance*leftWavelength) / (rightWavelength - leftWavelength);
continuum = gradient * spectrum.wavelength + intercept;
% Plot the continuum.
hold on;
plot(spectrum.wavelength, continuum, 'Color', 'r');
% Find the feature wavelength. This is the minimum value of the reflectance
% within the continuum interval.
% a is the start index of the continuum interval.
% b is the end index of the continuum interval.
intervalWavelength = spectrum.wavelength(spectrum.wavelength > leftWavelength & spectrum.wavelength < rightWavelength);
intervalReflectance = spectrum.data(spectrum.wavelength > leftWavelength & spectrum.wavelength < rightWavelength);
featureWavelength = intervalWavelength(intervalReflectance == min(intervalReflectance));
featureReflectance = spectrum.data(spectrum.wavelength == featureWavelength);
if numel(featureWavelength) > 1
error('Feature does not have a single minimum value.');
end
% Find the continuum reflectance. This is the reflectance at the
% feature wavelength.
continuumReflectance = gradient*featureWavelength + intercept;
% Plot a line to show the feature wavelength.
line([featureWavelength, featureWavelength], [featureReflectance, continuumReflectance], 'Color', 'c');
% Calculate the depth of the feature.
depth = 1 - (featureReflectance / continuumReflectance);
% Switch off graphics holding.
hold off;
end