# R/crcal.R In cryst: Calculate the Relative Crystallinity of Starch by XRD and FTIR

#### Documented in crystMWcrystSGfitFTIRcftirplotxrdplot

#' @title Relative Crystallinity Calculation of X-Ray Diffraction Pattern of Starch by Bruckner Method
#' @details Calculate the relative starch crystallinity of XRD pattern by Bruckner method.
#' @description Allow  to calculate the relative crystallinity of starch by XRD.
#' The basic concept of Bruckner approach involves obtaining a smoothed
#' line that separates the amorphous and crystalline starch of an X-ray
#' diffraction pattern. This smoothed line is achieved by applying
#' a moving average smoothing method to the original pattern.
#' @author Claudio Pozo Valenzuela [aut, cre] and Saddys Rodriguez-llamazares [aut]
#' @param pattern matrix. The matrix of X-ray diffraction pattern.  The first row corresponds to Bragg angle 2\eqn{\theta}; the second row corresponds to intensity.
#' @param N numeric. N length of the smoothing window (number of variables). Defaults to 41.
#' @param iter numeric. Iter number of iterations. Defaults to 100.
#' @return An object of class crystMW, which is a list with the following components:
#'   \item{original}{Original matrix of X-ray diffraction patterns.}
#'   \item{background}{Estimation of the background shape (curve of the amorphous starch).}
#'   \item{corrected}{Estimation of residual crystalline area (curve of the crystalline starch).}
#'   \item{summary}{Summary calculation of crystallinity. Total area under the curve of the diffraction pattern (A.U.); Amorphous area (A.U.); Crystalline area (A.U.); Relative crystallinity (\%).}
#' @importFrom flux auc
#' @examples
#' data(XRD)
#' # Convert data frame to matrix, select A-type starch
#' pattern <- as.matrix(t(XRD[,c("Bragg_angle","A")]))
#' # List of crystallinity components
#' crs <- crystMW(pattern, N = 41, iter = 100)
#' # Original matrix
#' original <- crs$original #' # Background shape #' background <- crs$background
#' # Curve of the crystalline starch
#' corrected <- crs$corrected #' # Summary calculation of crystallinity #' summary <- crs$summary
#' @references
#' Bruckner, S. (2000). "Estimation of the background in powder diffraction patterns through a robust smoothing procedure." Journal of Applied Crystallography 33(3 Part 2): 977-979.
#' @export
crystMW<-function(pattern, N=41, iter=100){
if (pattern[1, 1] < pattern[1, dim(pattern)[2]]) {
pattern <- t(apply(pattern, 1, rev))
}
if (missing(pattern)) {
stop('No pattern data provided')
}
if (missing(N)) {
stop('No length of the smoothing window provided')
}
if (missing(iter)) {
stop('No number of iterations provided')
}
# Inserting on both sides N points
s  <- pattern[2,]
sa <- rep(NA, length(s)+2*N)
sa <- c(rep(pattern[2,1],N),s,rep(pattern[2,length(s)],N))
# Iteration
for(it in 1:iter){
# Replacement of the kth point by the average of its 2N next neighbours (N to the left, N to the right)
Ik <- sa
for(k in N:(length(sa)-N)){
Ik[k] <- (1/(2*N))*(sum(sa[(k-N):(k-1)])+sum(sa[(k+1):(k+N)]))
}
# The new pattern is compared, point by point, with the original one, retaining the lowest of the two intensities
for(i in 1:length(sa)){
if(Ik[i] <= sa[i]){
sa[i] <- Ik[i]
}else{
sa[i] <- sa[i]
}
}
}
# Reconstructed background
background <-as.matrix(rbind(pattern[1,], sa[(N+1):(length(sa)-N)]))
colnames(background) <- NULL
rownames(background) <- NULL
# Corrected pattern (original minus background)
corrected<-as.matrix(rbind(pattern[1,], s-background[2,]))
colnames(corrected) <- NULL
rownames(corrected) <- NULL
# Crystalline area
res <- matrix(NA, nrow=1, ncol=4)
colnames(res) <-c('Total Area', 'Amorphous area', 'Crystalline area', 'Relative crystallinity')
rownames(res) <- NULL
## Total area
res[1] <- auc(pattern[1,], pattern[2,], thresh = 0, dens = 100)
## Amorphous area
res[2] <- auc(background[1,], background[2,], thresh = 0, dens = 100)
## Crystalline area
res[3] <-res[1]-res[2]
## Relative crystallinity
res[4] <- res[3]/res[1]
# Output
value <-list(original=pattern, background=background, corrected=corrected, summary=res)
attr(value, "class") <- "crystMW"
value
}

#' @title Relative Crystallinity Calculation of X-Ray Diffraction Pattern of Starch by Frost Method
#' @details Calculate the relative starch crystallinity of XRD pattern by Frost method.
#' @description Allow  to calculate the relative crystallinity of starch by XRD.
#' The basic concept of Frost approach involves obtaining a smoothed
#' line that separates the amorphous and crystalline starch of an X-ray
#' diffraction pattern. This smoothed line is achieved by applying
#' a Savitzky-Golay smoothing method to the original pattern.
#' @author Claudio Pozo Valenzuela [aut, cre] and Saddys Rodriguez-llamazares [aut]
#' @param pattern matrix. The matrix of X-ray diffraction pattern.  The first row corresponds to Bragg angle 2\eqn{\theta}; the second row corresponds to intensity.
#' @param N numeric. N length of the smoothing window (number of variables). Defaults to 101.
#' @param iter numeric. Iter number of iterations. Defaults to 400.
#' @param p numeric. Filter order. Defaults to 2.
#' @return An object of class crystSG, which is a list with the following components:
#'   \item{original}{Original matrix of X-ray diffraction patterns.}
#'   \item{background}{Estimation of the background shape (curve of the amorphous starch).}
#'   \item{corrected}{Estimation of residual crystalline area (curve of the crystalline starch).}
#'   \item{summary}{Summary calculation of crystallinity. Total area under the curve of the diffraction pattern (A.U.); Amorphous area (A.U.); Crystalline area (A.U.); Relative crystallinity (\%).}
#' @importFrom flux auc
#' @importFrom stats convolve
#' @importFrom pracma pinv
#' @examples
#' data(XRD)
#' # Convert data frame to matrix, select A-type starch
#' pattern <- as.matrix(t(XRD[, c("Bragg_angle","A")]))
#' # List of crystallinity components
#' crs <- crystSG(pattern, N = 101, iter = 400, p = 2)
#' # Original matrix
#' original <- crs$original #' # Background shape #' background <- crs$background
#' # Curve of the crystalline starch
#' corrected <- crs$corrected #' # Summary calculation of crystallinity #' summary <- crs$summary
#' @references
#' Frost, K., et al. (2009). "Crystallinity and structure of starch using wide angle X-ray scattering." Carbohydrate Polymers 78(3): 543-548.
#' @export
crystSG<-function(pattern, N=101, iter=400, p=2){
if (pattern[1, 1] < pattern[1, dim(pattern)[2]]) {
pattern <- t(apply(pattern, 1, rev))
}
if (missing(pattern)) {
stop('No pattern data provided')
}
if (missing(N)) {
stop('No length of the smoothing window provided')
}
if (missing(iter)) {
stop('No number of iterations provided')
}
if (missing(p)) {
stop('No filter order provided')
}
# Inserting on both sides N points
s  <- pattern[2,]
sa <- rep(NA, length(s)+2*N)
sa <- c(rep(pattern[2,1],N),s,rep(pattern[2,length(s)],N))
# Iteration
for(it in 1:iter){
# Replacement of the kth point by the savitzky golay smoothing of its 2N next neighbours (N to the left, N to the right)
fc <- (N - 1)/2
X  <- outer(-fc:fc, 0:p, FUN = "^")
Y  <- pinv(X)
t  <- convolve(sa, rev(Y[1,]), type = "o")
Ik <- t[(fc + 1):(length(t) - fc)]
# The new pattern is compared, point by point, with the original one, retaining the lowest of the two intensities
for(i in 1:length(sa)){
if(Ik[i]<=sa[i]){
sa[i]<-Ik[i]
}else{
sa[i]<-sa[i]
}
}
}
# Reconstructed background
background <-as.matrix(rbind(pattern[1,], sa[(N+1):(length(sa)-N)]))
colnames(background) <- NULL
rownames(pattern)<- NULL
# Corrected pattern (original minus background)
corrected<-as.matrix(rbind(pattern[1,], s-background[2,]))
colnames(corrected) <- NULL
rownames(corrected) <- NULL
# Crystalline area
res <- matrix(NA, nrow=1, ncol=4)
colnames(res) <-c('Total Area', 'Amorphous area', 'Crystalline area', 'Relative crystallinity')
rownames(res) <- NULL
## Total area
res[1] <- auc(pattern[1,], pattern[2,], thresh = 0, dens = 100)
## Amorphous area
res[2] <- auc(background[1,], background[2,], thresh = 0, dens = 100)
## Crystalline area
res[3] <-res[1]-res[2]
## Relative crystallinity
res[4] <- res[3]/res[1]
# Output
value <-list(original=pattern, background=background, corrected=corrected, summary=res)
attr(value, "class") <- "crystSG"
value
}

#' @title Relative Crystallinity Calculation of FTIR Spectrum of Starch by SUN Method
#' @details Calculate the relative starch crystallinity of FTIR spectrum by SUN method.
#' @description Allow  to calculate the relative crystallinity of starch by FTIR.
#' The basic concept of SUN approach involves obtaining a gaussian holocrystalline-peak
#' in the 800-1300 cm-1 region of FTIR spectrum of starch which is divided into amorphous
#' region and crystalline region.
#' @author Claudio Pozo Valenzuela [aut, cre] and Saddys Rodriguez-llamazares [aut]
#' @param spectrum matrix. The matrix of FTIR spectrum baseline-corrected by drawing a tangentline in the 800-1300 cm-1 region. The first row corresponds to wavelength; the second row corresponds to intensity.
#' @param mu numeric. Gaussian mean of holocrystalline-peak. Defaults to 1180.
#' @param sigma numeric. Standard deviation of holocrystalline-peak. Defaults to 60.
#' @param k numeric. Arbitrary scaling parameter. Defaults to 1.
#' @param lim vector. Fitting points of holocrystalline-peak. Defaults to c(1190, 1160, 985, 950).
#' @return An object of class fitFTIRc, which is a list with the following components:
#'   \item{original}{Original matrix of FTIR spectrum.}
#'   \item{gauss}{Gaussian curve fit.}
#'   \item{fit}{Summary of Non-Linear Least-Squares Model Fits.}
#'   \item{summary}{Summary calculation of crystallinity. Total area under the curve of the diffraction spectrum (A.U.); Amorphous area (A.U.); Crystalline area (A.U.); Relative crystallinity (\%).}
#' @importFrom flux auc
#' @importFrom stats nls
#' @importFrom stats coef
#' @examples
#' # Convert data frame to matrix, select A-type starch
#' spectrum <- as.matrix(t(FTIR[, c('wavelength','A')]))
#' # List of crystallinity components
#' crs <- fitFTIRc(spectrum = spectrum, mu = 1180, sigma = 60, k = 1, lim = c(1190, 1160, 985, 955))
#' # Original matrix
#' original <- crs$original #' # Gaussian curve fit #' gauss <- crs$gauss
#' # Summary of Non-Linear Least-Squares Model Fits
#' fit <- crs$fit #' # Summary calculation of crystallinity #' summary <- crs$summary
#' @references
#' Sun, Y., et al. (2014). "A new method for determining the relative crystallinity of chickpea starch by Fourier-transform infrared spectroscopy." Carbohydrate Polymers 108: 153-158.
#' @export
fitFTIRc<- function(spectrum, mu=1180, sigma=60, k=1, lim=c(1190, 1160, 985, 950)){
if (spectrum[1, 1] < spectrum[1, dim(spectrum)[2]]) {
spectrum <- t(apply(spectrum, 1, rev))
}
if (missing(spectrum)) {
stop('No spectrum data provided')
}
if (missing(mu)) {
stop('No gaussian mean provided')
}
if (missing(sigma)) {
stop('No standard deviation provided')
}
if (missing(k)) {
stop('No arbitrary scaling parameter provided')
}
if (missing(lim)) {
stop('No fitting points of holocrystalline-peak provided')
}
# Truncated region between 950-985 cm-1 and 1160-1190 cm-1
lim<-sort(lim, decreasing = TRUE)
li01 <- which(round(spectrum[1, ], 0) == lim[1])[1]
ls01 <- which(round(spectrum[1, ], 0) == lim[2])[1]
li02 <- which(round(spectrum[1, ], 0) == lim[3])[1]
ls02 <- which(round(spectrum[1, ], 0) == lim[4])[1]
a01 <- spectrum[, ls01:li01]
a02 <- spectrum[, ls02:li02]
a03 <- cbind(c(1300,rep(0, time=1)),a02,a01, c(800,rep(0, time=1)))
#Apply function nls
## First present the data in a data frame
s <- data.frame(x=a03[1,], r=a03[2,])
fit <- nls(r ~ (k*exp(-(x-mu)**2/(2*sigma**2))),data=s, start=list(k=2, mu=1063, sigma=50))
# Crystalline area
res <- matrix(NA, nrow=1, ncol=4)
colnames(res) <-c('Total Area', 'Amorphous area', 'Crystalline area', 'Relative crystallinity')
rownames(res) <- NULL
## Gaussian line polygon
gl <- matrix(NA, nrow=2, ncol=ncol(spectrum))
gl[1,] <- spectrum[1,]
for (i in 1:length(spectrum[1,])){
gl[2,i]<-coef(fit)[1] * exp(-(gl[1,i]-coef(fit)[2])**2/(2*coef(fit)[3]**2))
}
## Crystalline polygon
pol <- matrix(NA, nrow=2, ncol=ncol(spectrum))
pol[1,] <- spectrum[1,]
for (i in 1:length(spectrum[1,])){
if (gl[2,i]<=spectrum[2,i]){
pol[2,i]<-gl[2,i]
}else{
pol[2,i]<-spectrum[2,i]
}
}
## Total area
res[1] <- auc(gl[1,], gl[2,], thresh = 0, dens = 100) # Gaussian line polygon
## Crystalline area
res[2] <- auc(pol[1,], pol[2,], thresh = 0, dens = 100) # Crystalline polygon
## Amorphous area
res[3] <-res[1]-res[2]
## Relative crystallinity
res[4] <- res[2]/res[1]
# Output
value <- list(original=spectrum, gauss=gl, fit=summary(fit), summary=res)
attr(value, "class") <- "fitFTIRc"
value
}

#' @title Plots the Crystalline Area of a X-Ray Diffraction Pattern of Starch
#' @description Produces a graph of the crystalline area of a X-ray diffraction pattern of starch and background.
#' @author Claudio Pozo Valenzuela [aut, cre] and Saddys Rodriguez-llamazares [aut]
#' @param pattern matrix. The matrix of X-ray diffraction pattern.  The first row corresponds to Bragg angle 2\eqn{\theta}; the second row corresponds to intensity.
#' @param background matrix. The matrix of background shape (curve of the amorphous starch). The first row corresponds to Bragg angle 2\eqn{\theta}; the second row corresponds to intensity.
#' @importFrom graphics plot
#' @importFrom graphics lines
#' @importFrom graphics polygon
#' @importFrom graphics title
#' @importFrom grDevices rgb
#' @examples
#' # Convert data frame to matrix, select A-type starch
#' pattern <- as.matrix(t(XRD[, c("Bragg_angle","A")]))
#' # List of crystallinity components
#' crs <- crystMW(pattern, N = 11, iter = 100)
#' # Original matrix
#' original <- crs$original #' # Background shape #' background <- crs$background
#' # Plots the crystalline area of a XRD pattern
#' xrdplot(pattern=original, background=background)
#' @export
xrdplot <- function(pattern, background){
if (missing(pattern)) {
stop('No pattern data provided')
}
if (missing(background)) {
stop('No background provided')
}
if (pattern[1, 1] < pattern[1, dim(pattern)[2]]) {
pattern <- t(apply(pattern, 1, rev))
}
if (background[1, 1] < background[1, dim(background)[2]]) {
background <- t(apply(background, 1, rev))
}

titx= '2'~{theta}         # x title
tity= 'Intensity (a.u.)'  # y title
cexLab= 1.5

# Limits
xlimS <- round(max(pattern[1,]), 0)
xlimI <- round(min(pattern[1,]), 0)
ylimS <- round(max(pattern[2, ])*1.3,1)
ylimI <- round(min(background[2, ]),2)
xlim<-c(xlimI, xlimS)
ylim<-c(ylimI, ylimS)

# Plot XRD
plotpattern<-plot(pattern[1, ], pattern[2, ],  # function to plot
xlab="",          # suppress x labels
ylab="",          # suppress y labels
type = 'l',       # specify line graph
xlim = xlim,      # extend axis limits to give space for text annotation
ylim = ylim       # ditto
)
lines(background[1, ], background[2, ],col="black")
polygon(c(pattern[1,],rev(background[1,])),c(pattern[2,],rev(background[2,])),col="grey")

# Titles
title(xlab= titx, col.lab=rgb(0, 0, 0),  cex.lab=cexLab, line = 3.8)
title(ylab= tity, col.lab=rgb(0, 0, 0),  cex.lab=cexLab)
}

#' @title Plots the Crystalline Area of a FTIR Spectrum of Starch
#' @description Produces a graph of the crystalline area of a FTIR spectrum of starch and the Gauss curve.
#' @author Claudio Pozo Valenzuela [aut, cre] and Saddys Rodriguez-llamazares [aut]
#' @param spectrum matrix. The matrix of FTIR spectrum baseline-corrected by drawing a tangentline in the 800-1300 cm-1 region. The first row corresponds to wavelength; the second row corresponds to intensity.
#' @param gauss matrix. The matrix of Gauss curve (gaussian holocrystalline-peak).
#' @param lim vector. Regions of the FTIR spectrum comprising the fixing points of the Gauss curve.
#' @importFrom graphics plot
#' @importFrom graphics lines
#' @importFrom graphics abline
#' @importFrom graphics title
#' @importFrom graphics polygon
#' @importFrom grDevices rgb
#' @examples
#' # Convert data frame to matrix, select A-type starch
#' spectrum <- as.matrix(t(FTIR[, c('wavelength','A')]))
#' # List of crystallinity components
#' crs <- fitFTIRc(spectrum = spectrum, mu = 1180, sigma = 60, k = 1, lim = c(1190, 1160, 985, 955))
#' # Original matrix
#' original <- crs$original #' # Gaussian curve fit #' gauss <- crs$gauss
#' # Plots the crystalline area of a FTIR spectrum
#' ftirplot(spectrum=original, gauss=gauss, lim=c(1190, 1160, 985, 950))
#' @export
ftirplot <- function(spectrum, gauss, lim=c(1190, 1160, 985, 950)) {
if (missing(spectrum)) {
stop('No spectrum data provided')
}
if (missing(gauss)) {
stop('No gauss curve provided')
}
if (missing(lim)) {
stop('No fitting points of holocrystalline-peak provided')
}
if (spectrum[1, 1] < spectrum[1, dim(spectrum)[2]]) {
spectrum <- t(apply(spectrum, 1, rev))
}
if (gauss[1, 1] < gauss[1, dim(gauss)[2]]) {
gauss <- t(apply(gauss, 1, rev))
}
# Titles
titx= 'Wavelength'~(cm^{-1}) # x title
tity= "Absorbance (a.u.)"    # y title
cexLab= 1.5                  # magnification of x and y labels
# Limits
xlimS <- round(max(spectrum[1,]), -1)
xlimI <- round(min(spectrum[1,]), -1)
ylimS <- round(max(gauss[2,])*1.2,1)
ylimI <- round(min(spectrum[2, ]),2)
xlim<-c(xlimS, xlimI)
ylim<-c(ylimI, ylimS)
# Plot spectrum and abline
plot(gauss[1,], gauss[2,], type='l',xlab="", ylab="",yaxt='n', xlim=xlim, ylim=ylim)
lines(spectrum[1,],spectrum[2,],col="black")
lim<-sort(lim, decreasing = TRUE)
abline(v=lim[4], lty= 2)
abline(v=lim[3], lty= 2)
abline(v=lim[2], lty= 2)
abline(v=lim[1], lty= 2)
# Calculate polygon of crystalline starch
pol <- matrix(NA, nrow=2, ncol=ncol(spectrum))
pol[1,] <- spectrum[1,]
for (i in 1:length(spectrum[1,])){
if (gauss[2,i]<=spectrum[2,i]){
pol[2,i]<-gauss[2,i]
}else{
pol[2,i]<-spectrum[2,i]
}
}
# Crystalline polygon
pol <- matrix(NA, nrow=2, ncol=ncol(spectrum))
pol[1,] <- spectrum[1,]
for (i in 1:length(spectrum[1,])){
if (gauss[2,i]<=spectrum[2,i]){
pol[2,i]<-gauss[2,i]
}else{
pol[2,i]<-spectrum[2,i]
}
}
# Titles
title(xlab= titx, col.lab=rgb(0, 0, 0),  cex.lab=cexLab, line = 3.4)
title(ylab= tity, col.lab=rgb(0, 0, 0),  cex.lab=cexLab, line = 1)
# Apply cross-hatching to a polygon
polygon(pol[1, ], pol[2, ], density=10, angle=45, border=NULL)
}


## Try the cryst package in your browser

Any scripts or data that you put into this service are public.

cryst documentation built on May 29, 2017, 11 a.m.