#' Returns colors from colorbrewer2 palette
#'
#' @param ncols
#' number of colors
#'
#' @export
getColors = function(ncols = 1)
{
palette = c("#3288BD", "#66C2A5", "#ABDDA4", "#E6F598",
"#FEE08B", "#FDAE61", "#F46D43", "#D53E4F")
if (ncols == 3)
{
colvals = c("#3288BD", "#91CF60", "#D53E4F")
}
else
{
colfunc = colorRampPalette(palette, space = 'Lab')
colvals = colfunc(ncols)
}
colvals
}
#' Normalization of spectra
#'
#' @description
#' normalize spectra to unit area
#'
#' @param spectra
#' matrix of vector with spectral data
#'
#' @export
areanorm = function(spectra)
{
spectra = as.matrix(spectra)
w = rowSums(abs(spectra))
normspectra = sweep(spectra, 1, w, '/')
normspectra
}
#' Purity spectra plot
#'
#' @description
#' Shows a plot with purity spectra for each or user selected component
#'
#' @param obj
#' object with unmixing results
#' @param ncomp
#' which component to make the plot for
#' @param origspec
#' original spectral data used for unmixing
#' @param col
#' colors for the std, purity and original spectra
#' @param lty
#' line type for the component spectra
#' @param lwd
#' line width for the component spectra
#' @param xlab
#' label for x axis
#' @param ylab
#' label for y axis
#' @param main
#' main title for the plot
#'
#' @export
plotPurity = function(obj, ncomp = 1, origspec = NULL,
col = c('#62B7EF', '#D84E5E', '#F0F0F0'), lty = 1, lwd = 1,
xlab = 'Wavenumbers', ylab = c('Weighted std', 'Purity'),
main = sprintf('Purity for C%d', ncomp))
{
if (class(obj) != 'purity')
stop('This plot can be maid only for "purity" objects!')
if (length(col) != 3)
stop('Parameter "col" should have three values!')
if (length(ncomp) != 1)
stop('Parameter "ncomp" should have one value!')
purityspec = obj$purityspec[ncomp, , drop = F]
stdspec = obj$stdspec[ncomp, , drop = F]
if (!obj$by.spec)
wavelength = 1:nrow(obj$conc)
else
wavelength = obj$wavelength
xlim = c(min(wavelength), max(wavelength))
ylim = c(min(stdspec), max(stdspec))
dx = xlim[2] - xlim[1]
dy = ylim[2] - ylim[1]
# correct margin to have two y-axis
par(mar = c(5, 4, 4, 5) + .1)
# first plot with std and spectra
plot(0,
type = 'n',
ylab = ylab[1],
xlab = xlab,
main = main,
xlim = c(xlim[1], xlim[2] + dx/20),
ylim = c(ylim[1], ylim[2] + dy/20)
)
# if original data is provided scale and show it
if (!is.null(origspec))
{
if (!obj$by.spec)
origspec = t(origspec)
origspec = (origspec - min(origspec))/(max(origspec) - min(origspec))
origspec = ylim[1] + origspec * (ylim[2] - ylim[1])
matlines(wavelength, t(origspec), type = 'l', col = col[3], lty = lty, lwd = lwd)
}
# show standard deviation plot
lines(wavelength, t(stdspec), lty = lty, lwd = lwd, col = col[1])
# show legend
l = legend(0, 0, c('Std', 'Purity'), col = col, cex = 0.80, lty = lty, lwd = lwd, plot = F)
legend(xlim[2] - l$rect$w + dx/20, ylim[2] + dy/20, c('Std', 'Purity'), col = col,
lty = lty, lwd = lwd, cex = 0.80, box.col = '#909090', bg = '#ffffff')
# second plot with purity spectrum
par(new = TRUE)
matplot(wavelength, t(purityspec), type = 'l', lty = lty, lwd = lwd, col = col[2],
xlim = c(xlim[1], xlim[2] + dx/20),
xaxt = 'n',
yaxt = 'n', xlab = '', ylab = '')
axis(4)
mtext(ylab[2], side = 4, line = 3, cex = 0.7)
points(wavelength[obj$purevars[ncomp]], purityspec[obj$purevars[ncomp]], pch = 16, col = col[2])
}
#' Spectral plot
#'
#' @description
#' Shows a plot with resolved spectra
#'
#' @param obj
#' object with unmixing results
#' @param comp
#' which components make the plot for
#' @param origspec
#' original spectral data used for unmixing
#' @param show.labels
#' logical, show or not labels for the selected pure variables on the plot
#' @param show.legend
#' logical, show or not a legend for the plot
#' @param legend.str
#' a vector with legend items
#' @param xlim
#' limits for x-axis
#' @param ylim
#' limits for y-axis
#' @param col
#' color for the component spectra
#' @param lty
#' line type for the component spectra
#' @param lwd
#' line width for the component spectra
#' @param xlab
#' label for x axis
#' @param ylab
#' label for y axis
#' @param main
#' main title for the plot
#' @param ...
#' other maplot parameters
#'
#' @export
plotSpectra = function(obj, comp = 1:obj$ncomp, origspec = NULL, show.labels = T,
show.legend = T, legend.str = NULL, ylim = NULL, xlim = NULL,
col = getColors(length(comp)), lty = 1, lwd = 1,
xlab = 'Wavenumbers', ylab = 'Intensity',
main = 'Resolved spectra', ...)
{
purespec = areanorm(t(obj$spec[, comp, drop = F]))
purevars = obj$purevars
wavelength = obj$wavelength
if (!is.null(origspec))
{
origspec = areanorm(origspec)
data = rbind(origspec, purespec)
}
else
{
data = rbind(purespec)
}
if (is.null(xlim))
xlim = c(min(wavelength), max(wavelength))
if (is.null(ylim))
ylim = c(min(data), max(data))
dx = xlim[2] - xlim[1]
dy = ylim[2] - ylim[1]
if (!is.null(origspec))
matplot(wavelength, t(origspec), type = 'l', col = '#F0F0F0', lty = lty, lwd = lwd,
ylab = ylab,
xlab = xlab,
main = main,
xlim = c(xlim[1] - dx/20, xlim[2] ),
ylim = c(ylim[1], ylim[2] + dy/20),
...)
else
matplot(0, 0, type = 'n', col = '#F0F0F0', lty = lty, lwd = lwd,
ylab = ylab,
xlab = xlab,
main = main,
xlim = c(xlim[1] - dx/20, xlim[2] ),
ylim = c(ylim[1], ylim[2] + dy/20),
...)
matlines(wavelength, t(purespec), lty = lty, lwd = lwd, col = col)
if (show.legend)
{
if (is.null(legend.str))
legend.str = colnames(obj$spec)[comp]
legend(xlim[1] - dx/16, ylim[2] + dy/25, legend.str, col = col, lty = lty,
lwd = lwd, cex = 0.85, box.col = '#909090')
}
if (show.labels & !is.null(purevars))
{
for (i in 1:length(comp))
{
text(wavelength[purevars[comp[i]]], purespec[i, purevars[comp[i]]],
sprintf('p%d', comp[i]), cex = 0.85, pos = 3)
}
}
}
#' Concentrations plot
#'
#' @description
#' Shows a plot with resolved concentration profiles
#'
#' @param obj
#' object with unmixing results
#' @param comp
#' which components make the plot for
#' @param origspec
#' original spectral data used for unmixing
#' @param col
#' color for the plot lines
#' @param lty
#' line type for plot lines
#' @param lwd
#' line width for plot lines
#' @param xlab
#' label for x axis
#' @param ylab
#' label for y axis
#' @param main
#' main title for the plot
#' @param ...
#' other parameters for matplot function
#'
#' @export
plotConcentrations = function(obj, comp = 1:obj$ncomp, origspec = NULL,
col = getColors(length(comp)), lty = 1, lwd = 1,
xlab = 'Mixings', ylab = '',
main = 'Resolved concentrations', ...)
{
pureconc = areanorm(t(obj$conc[, comp, drop = F]))
x = 1:nrow(obj$conc)
if (!is.null(origspec))
{
origconc = areanorm(t(origspec))
data = rbind(origconc, pureconc)
}
else
{
origconc = NULL
data = pureconc
}
xlim = c(min(x), max(x))
ylim = c(min(data), max(data))
dx = xlim[2] - xlim[1]
dy = ylim[2] - ylim[1]
if (!is.null(origconc))
matplot(x, t(origconc), type = 'l', col = '#F0F0F0', lty = lty, lwd = lwd,
ylab = ylab,
xlab = xlab,
main = main,
xlim = c(xlim[1] - dx/20, xlim[2] ),
ylim = c(ylim[1], ylim[2] + dy/20),
...)
else
matplot(0, 0, type = 'n', col = '#F0F0F0', lty = lty, lwd = lwd,
ylab = ylab,
xlab = xlab,
main = main,
xlim = c(xlim[1] - dx/20, xlim[2] ),
ylim = c(ylim[1], ylim[2] + dy/20),
...)
matlines(x, t(pureconc), lty = lty, lwd = lwd, col = col)
legend(xlim[1] - dx/16, ylim[2] + dy/25, colnames(obj$spec)[comp], col = col, lty = lty,
lwd = lwd, cex = 0.85, box.col = '#909090')
}
#' Residual variance plot for results of unmixing
#'
#' @param obj
#' object with unmixing results
#' @param type
#' type of the plot
#' @param col
#' color for line and markers
#' @param pch
#' marker symbol
#' @param xlab
#' label for x axis
#' @param ylab
#' label for y axis
#' @param main
#' main title for the plot
#' @param ...
#' other parameters for plot function
#'
#' @export
plotVariance = function(obj, type = 'b', col = getColors(1), pch = 16,
xlab = 'Number of components', ylab = 'Variance, %',
main = 'Residual variance', ...)
{
plot(1:obj$ncomp, 100 * obj$resvar, type = type, col = col, pch = pch, xlab = xlab, ylab = ylab,
main = main, xaxt = 'n', ...)
axis(1, at = 1:obj$ncomp, labels = 1:obj$ncomp)
}
#' Compares resolved and original pure component spectra
#'
#' @description
#' shows two spectra - a resolved and an original for selected pure component in the same plot
#'
#' @param obj
#' object with unmixing results
#' @param purespec
#' matrix with pure spectra for the components (spectra in rows)
#' @param ncomp
#' for which component show the comparison for
#' @param col
#' color for the spectra (two values - one for each spectrum)
#' @param lwd
#' line width for the spectra (two values)
#' @param lty
#' line type for the spectra (two values)
#' @param xlab
#' label for x axis
#' @param ylab
#' label for y axis
#' @param main
#' main title for the plot
#' @param show.stat
#' logical, show or not relative residual variance between the spectra
#' @param ...
#' other plot parameters
#' @export
plotSpectraComparison = function(obj, purespec, ncomp = 1, col = getColors(2), lwd = c(2, 1),
lty = c(1, 1), xlab = 'Wavenumbers', ylab = 'Intensity',
main = NULL, show.stat = F, ...)
{
rspectrum = t(obj$spec[, ncomp, drop = F])
ospectrum = t(purespec[, ncomp, drop = F])
ospectrum = areanorm(ospectrum)
rspectrum = areanorm(rspectrum)
res = sum((ospectrum - rspectrum)^2) / sum(ospectrum^2)
if (is.null(main))
{
if (show.stat)
main = sprintf('Spectra comparison for C%d (%.3f%%)', ncomp, 100 * res)
else
main = sprintf('Spectra comparison for C%d', ncomp)
}
xlim = c(min(obj$wavelength), max(obj$wavelength))
ylim = c(min(c(ospectrum, rspectrum)), max(c(ospectrum, rspectrum)))
dx = xlim[2] - xlim[1]
dy = ylim[2] - ylim[1]
plot(obj$wavelength, ospectrum, type = 'l', main = main, xlab = xlab, ylab = ylab,
xlim = c(xlim[1] - dx/20, xlim[2]), ylim = c(ylim[1], ylim[2] + dy/20),
col = col[1], lty = lty[1], lwd = lwd[1], ...)
lines(obj$wavelength, rspectrum, col = col[2], lty = lty[2], lwd = lwd[2])
abline(v = obj$wavelength[obj$purevars[ncomp]], lty = 2, lwd = 2, col = '#909090')
legend(xlim[1] - dx/15, ylim[2] + dy/25, c('Original', 'Resolved'), col = col, lty = lty,
lwd = lwd, cex = 0.80, box.col = '#909090')
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.