Skip to content

Commit

Permalink
- Improved smoothing and plot after processing.
Browse files Browse the repository at this point in the history
  • Loading branch information
@LlucSF
LlucSF committed Nov 19, 2019
1 parent 4be5348 commit df12ab5
Show file tree
Hide file tree
Showing 2 changed files with 60 additions and 59 deletions.
43 changes: 10 additions & 33 deletions R/AverageSpectrum.R
View file
Original file line number Diff line number Diff line change
Expand Up @@ -16,60 +16,37 @@
# along with this program. If not, see <http://www.gnu.org/licenses/>.
############################################################################

AverageSpectrum <- function(numBands, numPixels, rCRSIObj)
AverageSpectrum <- function(numBands, numPixels, rCRSIObj, plot = T)
{
temp<- array(dim =(c(rCRSIObj$BandsLength,numBands)))

for (i in 1:numBands)
{
for (j in 1:rCRSIObj$BandsLength)
{
temp[j,i] = 0
temp[j,i] <- 0
}
}

for (i in 1:numBands)
{
for (j in 1:numPixels)
{
temp[,i] = temp[,i] + rCRSIObj$Data[(j+(numPixels*(i-1))),]
temp[,i] <- temp[,i] + rCRSIObj$Data[(j+(numPixels*(i-1))),]
}
temp[,i] = temp[,i] / numPixels
temp[,i] <- temp[,i] / numPixels
}

for (i in 1:numBands)
{
g <- ggplot2::ggplot(mapping = ggplot2::aes(y = temp[,i], x = rCRSIObj$RamanShiftAxis[,i])) + ggplot2::geom_line(colour = "red") +
ggplot2::ylim(c(min(temp[,i]), max(temp[,i]))) + ggplot2::ggtitle(label = "Average Spectrum") + ggplot2::labs(x ="Raman Shift", y = "Counts") +
ggplot2::theme_bw() + ggplot2::scale_x_continuous(breaks = trunc(seq(from = min(rCRSIObj$RamanShiftAxis[,i]), to = max(rCRSIObj$RamanShiftAxis[,i]),length.out = 30)),minor_breaks = ggplot2::waiver())
print(g)
}

for (i in 1:numBands)
if(plot)
{
g <- ggplot2::ggplot(mapping = ggplot2::aes(y = apply(rCRSIObj$Data[(numPixels*(i-1)+1):(numPixels*(i)),],2,max), x = rCRSIObj$RamanShiftAxis[,i])) + ggplot2::geom_line(colour = "red") +
ggplot2::ylim(c(min(apply(rCRSIObj$Data[(numPixels*(i-1)+1):(numPixels*(i)),],2,max)), max(apply(rCRSIObj$Data[(numPixels*(i-1)+1):(numPixels*(i)),],2,max)))) + ggplot2::ggtitle(label = "Skyline Spectrum") + ggplot2::labs(x ="Raman Shift", y = "Counts") +
ggplot2::theme_bw() + ggplot2::scale_x_continuous(breaks = trunc(seq(from = min(rCRSIObj$RamanShiftAxis[,i]), to = max(rCRSIObj$RamanShiftAxis[,i]),length.out = 30)),minor_breaks = ggplot2::waiver())
print(g)
}

if(!is.null(rCRSIObj$ProAvrgSpectr))
for (i in 1:numBands)
{
for (i in 1:numBands)
{
df <- data.frame(y = c(rCRSIObj$AvrgSpectr[,i]/max(rCRSIObj$AvrgSpectr[,i]), temp[,i]/max(temp[,i])),
x = rep(rCRSIObj$RamanShiftAxis[,i], times = 2),
cl = c(rep("Raw", times = length(temp[,i]) ),rep("Processed", times = length(temp[,i]))))

g <- ggplot2::ggplot() + ggplot2::theme_bw() +
ggplot2::geom_line(mapping = ggplot2::aes(x = df$x, y = df$y, colour = df$cl)) +
ggplot2::labs(x ="Raman Shift", y = "Counts") +
ggplot2::ggtitle("Data processing normalized results") +
ggplot2::labs(colour = "Data") +
ggplot2::scale_x_continuous(breaks = trunc(seq(from = min(df$x), to = max(df$x),length.out = 30)),minor_breaks = ggplot2::waiver())
g <- ggplot2::ggplot(mapping = ggplot2::aes(y = temp[,i], x = rCRSIObj$RamanShiftAxis[,i])) + ggplot2::geom_line(colour = "red") +
ggplot2::ylim(c(min(temp[,i]), max(temp[,i]))) + ggplot2::ggtitle(label = "Average Spectrum") + ggplot2::labs(x ="Raman Shift", y = "Counts") +
ggplot2::theme_bw() + ggplot2::scale_x_continuous(breaks = trunc(seq(from = min(rCRSIObj$RamanShiftAxis[,i]), to = max(rCRSIObj$RamanShiftAxis[,i]),length.out = 30)),minor_breaks = ggplot2::waiver())
print(g)
}
}
}

return(temp)
}
Expand Down
76 changes: 50 additions & 26 deletions R/Smoothing.R
View file
Original file line number Diff line number Diff line change
Expand Up @@ -35,14 +35,14 @@ SmoothSpectra <- function(rCRSIObj,method = "normal", OddfiltLength = 7,
BackgroundSubs = TRUE, window = 32
)
{

if(method!="normal" & method!="wavelets")
{
stop("Method must be 'normal' or 'wavelets'")
}

writeLines("Starting processing...")
start_time <- Sys.time()
OldrCRSIObj <- rCRSIObj

if(method=="wavelets")
{
Expand All @@ -56,12 +56,12 @@ SmoothSpectra <- function(rCRSIObj,method = "normal", OddfiltLength = 7,

Wav_c <- wavelets::modwt(X = spectrum, n.levels = 4, filter = "d4",
boundary = "periodic", fast = TRUE)
W<-Wav_c@W
W$W1<-W$W1-W$W1
W$W2<-W$W2-W$W2
W$W3<-W$W3-W$W3
W$W4<-W$W4-W$W4
Wav_c@W<-W
W <- Wav_c@W
W$W1 <- W$W1-W$W1
W$W2 <- W$W2-W$W2
W$W3 <- W$W3-W$W3
# W$W4 <- W$W4-W$W4
Wav_c@W <- W

rCRSIObj$Data[b+(a-1)*rCRSIObj$numPixels,]<-wavelets::imodwt(Wav_c)[(window+1):(window+rCRSIObj$BandsLength)]

Expand All @@ -71,21 +71,31 @@ SmoothSpectra <- function(rCRSIObj,method = "normal", OddfiltLength = 7,
rCRSIObj$Data[(b+((a-1)*rCRSIObj$numPixels)),],
rep(x = rCRSIObj$Data[(b+((a-1)*rCRSIObj$numPixels)),rCRSIObj$BandsLength],times = window))

Wav_c <- wavelets::modwt(X = spectrum, n.levels=7, filter="d6",
Wav_c <- wavelets::modwt(X = spectrum, n.levels=9, filter="d6",
boundary="periodic", fast=TRUE)
W<-Wav_c@W
W$W7<-W$W7-W$W7
Wav_c@W<-W
V<-Wav_c@V

V$V1[,1]<-V$V1[,1]-V$V1[,1]
V$V2[,1]<-V$V2[,1]-V$V2[,1]
V$V3[,1]<-V$V3[,1]-V$V3[,1]
V$V4[,1]<-V$V4[,1]-V$V4[,1]
V$V5[,1]<-V$V5[,1]-V$V5[,1]
V$V6[,1]<-V$V6[,1]-V$V6[,1]
V$V7[,1]<-V$V7[,1]-V$V7[,1]
Wav_c@V<-V
W <- Wav_c@W
# W$W1 <- W$W1-W$W1
# W$W2 <- W$W2-W$W2
# W$W3 <- W$W3-W$W3
W$W4 <- W$W4-W$W4
W$W5 <- W$W5-W$W5
W$W6 <- W$W6-W$W6
W$W7 <- W$W7-W$W7
W$W8 <- W$W8-W$W8
W$W9 <- W$W9-W$W9
Wav_c@W <- W

V <- Wav_c@V
V$V1[,1] <- V$V1[,1]-V$V1[,1]
V$V2[,1] <- V$V2[,1]-V$V2[,1]
V$V3[,1] <- V$V3[,1]-V$V3[,1]
V$V4[,1] <- V$V4[,1]-V$V4[,1]
V$V5[,1] <- V$V5[,1]-V$V5[,1]
V$V6[,1] <- V$V6[,1]-V$V6[,1]
V$V7[,1] <- V$V7[,1]-V$V7[,1]
V$V8[,1] <- V$V8[,1]-V$V8[,1]
V$V9[,1] <- V$V9[,1]-V$V9[,1]
Wav_c@V <- V

rCRSIObj$Data[b+(a-1)*rCRSIObj$numPixels,] <- wavelets::imodwt(Wav_c)[(window+1):(window+rCRSIObj$BandsLength)]
}
Expand All @@ -108,20 +118,34 @@ SmoothSpectra <- function(rCRSIObj,method = "normal", OddfiltLength = 7,
}
}
}

rCRSIObj$ProAvrgSpectr <- 1
rCRSIObj$ProAvrgSpectr <- AverageSpectrum(rCRSIObj$numBands,rCRSIObj$numPixels,rCRSIObj)
}


if(method=="normal")
{
for(i in 1:rCRSIObj$numSpectr)
{
rCRSIObj$Data[i,] <- signal::sgolayfilt(x = rCRSIObj$Data[i,],p = 2,n = OddfiltLength)
}
}

rCRSIObj$ProAvrgSpectr <- AverageSpectrum(rCRSIObj$numBands,rCRSIObj$numPixels,rCRSIObj,F)
Old_average <- AverageSpectrum(OldrCRSIObj$numBands, OldrCRSIObj$numPixels, OldrCRSIObj, F)

rCRSIObj$ProAvrgSpectr <- 1
rCRSIObj$ProAvrgSpectr <- AverageSpectrum(rCRSIObj$numBands,rCRSIObj$numPixels,rCRSIObj)
for (i in 1:rCRSIObj$numBands)
{
df <- data.frame(y = c(Old_average[,i]/max(Old_average[,i]),
rCRSIObj$ProAvrgSpectr[,i]/max(rCRSIObj$ProAvrgSpectr[,i])),
x = rep(rCRSIObj$RamanShiftAxis[,i], times = 2),
cl = rep(c("Raw","Processed"), each = length(rCRSIObj$ProAvrgSpectr[,i])))

g <- ggplot2::ggplot() + ggplot2::theme_bw() +
ggplot2::geom_line(mapping = ggplot2::aes(x = df$x, y = df$y, colour = df$cl)) +
ggplot2::labs(x ="Raman Shift", y = "Counts") +
ggplot2::ggtitle("Data processing normalized results") +
ggplot2::labs(colour = "Data") +
ggplot2::scale_x_continuous(breaks = trunc(seq(from = min(df$x), to = max(df$x),length.out = 30)),minor_breaks = ggplot2::waiver())
print(g)
}

end_time <- Sys.time()
Expand Down

0 comments on commit df12ab5

Please sign in to comment.