Options to calculate Reflectance

• Agustin.Lobo@geo3bcn.csic.es
• 20220925

This script compares 4 options to calculate reflectance from radiance images acquired by hyperspectral scanning systems, in particular Specim FX17 and FX10 mounted on desktop-size Lumo Scanner stands:

• Option 1. Radiance image divided by the global average spectra of the radiance white image
  Refl(i,j,k) = Rad(i,j,k)/Rad(k)
  
• Option 2. As (1), but subtracting the average spectra of the Dark image
  Refl(i,j,k) = (Rad(i,j,k)-Dark(k))/(Rad(k)-Dark(k))
  
• Option 3. Radiance image divided by the radiance value of the corresponding column in the white image
  Refl(i,j,k) = Rad(i,j,k)/Rad(j,k)
  
• Option 4. As (3), but subtracting the average spectra of the Dark image
  Refl(i,j,k) = (Rad(i,j,k)-Dark(k))/(Rad(j,k)-Dark(k))

Options 3 and 4 account for heterogeneous illumination along columns of the captured row, while options 2 and 4 account for dark noise.

#Directories
defaultW <- getOption("warn")
options(warn=-1)
RSpectDir <- "/home/alobo/owncloudRSpect/RSpect"
imadir <- "../FX17_DeceptionIsland31_2018-04-20_07-46-53/FX17/capture"
vdir <- "../FX17_DeceptionIsland31_2018-04-20_07-46-53/QGIS"
#R packages
library(terra, quietly=TRUE)
library(tidyterra, quietly=TRUE)
library(reshape2, quietly=TRUE)
library(ggplot2, quietly=TRUE)
#Required information
load(file.path(RSpectDir,"FX17bands.rda")) #central wavelength of FX17 spectral bands 
load(file.path(RSpectDir,"FX17weirds.rda"))#systematic abnormal values of our FX17 instrument
load(file.path(RSpectDir,"RefMCAA01.rda")) #certified reflectance spectra of WCS-MC-010 

Introduction

Keep in mind that these images are in fact stacks of 1 row x 650 column hyperspectral pixels, and that the lighting system does not provide a uniform illumination. The Radiance White image corresponds to a narrow white stick. An example:

  options(warn=-1)
  terraOptions(progress=0)
  fnameimaWref    <- list.files(imadir,patt=glob2rx("WHITEREF*.raw"))[1] 
  fnameimadarkref <- list.files(imadir,patt=glob2rx("DARKREF*.raw"))[1]
  fnameimaRad     <- list.files(imadir,patt=glob2rx("FX17*.raw"))[1]
  WIma <- rast(file.path(imadir,fnameimaWref))
  names(WIma) <- paste0("B",1:nlyr(WIma))
  suf <- unlist(strsplit(fnameimaWref,"WHITEREF_FX17_"))[2]
  suf <- unlist(strsplit(suf,".raw"))[1]
  fnameimaWref2   <- paste0("WHITEREF_BIG_FX17_",suf)
  WIma2 <- rast(file.path(imadir,fnameimaWref2))
  names(WIma2) <- paste0("B",1:nlyr(WIma))
  defbands=c(170, 14, 86)

  a <- subset(WIma2, subset=defbands)
  delme <- crop(a,ext(WIma))
  as1 <- stretch(delme,smin=0, smax=4096)
  as2 <- stretch(delme,minq=0.02, maxq=0.98)
  crs(as2) <- ""
  ggW2st <- ggplot() + 
    geom_spatraster_rgb(data = as2) +
    coord_fixed(ratio = 1) +
    theme_void() +  theme(plot.title = element_text(hjust = 0.125), plot.margin=unit(c(0,0,0,0),"mm")) +
    ggtitle(paste("Radiance White image \n(Bands ", paste(defbands,collapse=", "), "with 2% stretching)"))
  rm(a, as1, as2,delme)
  print(ggW2st)

options(warn=-1)
WIma2.df  <- as.data.frame(mean(WIma2[1,,,drop=FALSE]),  cell=TRUE,na.rm=FALSE)
meamWIma2.df <- data.frame(Column=WIma2.df$cell,
                             Radiance=c(WIma2.df$mean))
ggW1 <- ggplot(data=meamWIma2.df) +
    geom_point( aes(x=Column, y=Radiance),size=0.5) +
    xlab("Column") + ylab("Mean Band Radiance") +
    theme(legend.position = c(0.75, 0.15)) + theme(plot.title = element_text(hjust = 0.5)) +
    scale_size_identity() 
    #ggtitle("Radiance of White Reference")
  print(ggW1)

Reflectance of White image

Option 1

Radiance image divided by the global average spectra of the radiance white image

options(warn=-1)
ReflWIma1 <- WIma2/unlist(global(WIma2,mean,na.rm=TRUE))

Option 2

Radiance image divided by the global average spectra of the radiance white image, and subtracting the average spectra of the Dark image

options(warn=-1)
DarkIma   <- rast(file.path(imadir,fnameimadarkref))
d <- global(DarkIma, mean, na.rm=TRUE)
ReflWIma2 <- (WIma2-unlist(d))/(unlist(global(WIma2,mean,na.rm=TRUE))-unlist(d))

Option 3

Radiance image divided by the radiance value of the corresponding column in the white image

options(warn=-1)
ReflWIma3 <- WIma2/WIma2

Option 4

Radiance image divided by the radiance value of the corresponding column in the white image, with noise subtracted. (Obviously, for the case of the White image, options 3 and 4 are identical).

options(warn=-1)
ReflWIma4 <- (WIma2-unlist(d))/(WIma2-unlist(d))

Comparison

options(warn=-1)
ReflWIma1.df  <- as.data.frame(mean(ReflWIma1[1,,,drop=FALSE]),  cell=TRUE,na.rm=FALSE)
ReflWIma2.df  <- as.data.frame(mean(ReflWIma2[1,,,drop=FALSE]),  cell=TRUE,na.rm=FALSE)
ReflWIma3.df  <- as.data.frame(mean(ReflWIma3[1,,,drop=FALSE]),  cell=TRUE,na.rm=FALSE)
ReflWIma4.df  <- as.data.frame(mean(ReflWIma4[1,,,drop=FALSE]),  cell=TRUE,na.rm=FALSE)
meamReflWImas.df <- data.frame(Column=ReflWIma1.df$cell,
                             Source=rep(paste0("Option ", 1:4), rep(ncol(ReflWIma1), 4)),
                             Reflectance=c(ReflWIma1.df$mean,ReflWIma2.df$mean,
                                           ReflWIma3.df$mean,ReflWIma4.df$mean))
ggW1_W4 <- ggplot(data=meamReflWImas.df) +
    geom_point( aes(x=Column, y=Reflectance, color=Source),size=0.5) +
    xlab("Column") + ylab("Mean Reflectance across bands") +
    theme(legend.position = c(0.75, 0.15)) + theme(plot.title = element_text(hjust = 0.5)) +
    scale_size_identity() +
    ggtitle("Reflectance of White Reference")
  print(ggW1_W4)

Note that, in this rather extreme case, the illumination heterogeneity will likely imply a severe error in the calculated reflectance image if lighting differences along the white reference bar are not taken into account (options 1 and 2).

Reflectance of standard reference

In this section we use the radiance image of a Labsphere diffuse reflectance standard (WCS-MC-020), and show the reflectance spectra following all four options compared to the certified spectrum. We crop the object from the same scanned radiance dataset as above, for which we use an interactively digitized vector object.

options(warn=-1)
fnameimaRad <- list.files(imadir,patt=glob2rx("FX17*.raw"))[1]
RadIma <- rast(file.path(imadir,fnameimaRad))
names(RadIma) <- paste0("B",1:nlyr(RadIma))
v <- vect(file.path(vdir,"DI31.shp"))[6,]
ext(RadIma)[3:4] <- c(-ext(RadIma)[4],ext(RadIma)[3])#adapt to v geometry

a <- subset(RadIma, subset=defbands)
delme <- crop(a,ext(RadIma))
as1 <- stretch(delme,smin=0, smax=4096)
as2 <- stretch(delme,minq=0.02, maxq=0.98)
plotRGB(as2, ext=ext(v) + c(10,10,10,10) )
plot(v,add=TRUE)

rm(a, delme, as1, as2)

We calculate and plot Reflectance according to all four options and compare them to the certified spectrum:

options(warn=-1)
ext(RadIma)[3:4] <- c(ext(RadIma)[4],-ext(RadIma)[3])#recover raster geometry
ReflIma1 <- RadIma/unlist(global(WIma2,mean,na.rm=TRUE))
#DarkIma   <- rast(file.path(imadir,fnameimadarkref))
#d <- global(DarkIma, mean, na.rm=TRUE)
ReflIma2 <- (RadIma-unlist(d))/(unlist(global(WIma2,mean,na.rm=TRUE))-unlist(d))
ReflIma3 <- RadIma/WIma2
ReflIma4 <- (RadIma-unlist(d))/(WIma2-unlist(d))
#Extract values and calculate
ext(ReflIma1)[3:4] <- c(-ext(ReflIma1)[4],ext(ReflIma1)[3])#adapt to v geometry
ext(ReflIma2)[3:4] <- c(-ext(ReflIma2)[4],ext(ReflIma2)[3])#adapt to v geometry
ext(ReflIma3)[3:4] <- c(-ext(ReflIma3)[4],ext(ReflIma3)[3])#adapt to v geometry
ext(ReflIma4)[3:4] <- c(-ext(ReflIma4)[4],ext(ReflIma4)[3])#adapt to v geometry

refl1dat.m  <- terra::extract(ReflIma1, v, fun=mean, na.rm=TRUE)
refl1dat.sd <- terra::extract(ReflIma1, v, fun=sd, na.rm=TRUE)
refl2dat.m  <- terra::extract(ReflIma2, v, fun=mean, na.rm=TRUE)
refl2dat.sd <- terra::extract(ReflIma2, v, fun=sd, na.rm=TRUE)
refl3dat.m  <- terra::extract(ReflIma3, v, fun=mean, na.rm=TRUE)
refl3dat.sd <- terra::extract(ReflIma3, v, fun=sd, na.rm=TRUE)
refl4dat.m  <- terra::extract(ReflIma4, v, fun=mean, na.rm=TRUE)
refl4dat.sd <- terra::extract(ReflIma4, v, fun=sd, na.rm=TRUE)

refl1dat <- data.frame(Reflectance=t(refl1dat.m[-1]),sd=t(refl1dat.sd[-1]))
refl1dat$Bandnb <- FX17bands$Bandnb
refl1dat$Wavelength <- FX17bands$Wavelength
refl1dat$source <- "Option 1"
#head(refl1dat)

refl2dat <- data.frame(Reflectance=t(refl2dat.m[-1]),sd=t(refl2dat.sd[-1]))
refl2dat$Bandnb <- FX17bands$Bandnb
refl2dat$Wavelength <- FX17bands$Wavelength
refl2dat$source <- "Option 2"
#head(refl2dat)

refl3dat <- data.frame(Reflectance=t(refl3dat.m[-1]),sd=t(refl3dat.sd[-1]))
refl3dat$Bandnb <- FX17bands$Bandnb
refl3dat$Wavelength <- FX17bands$Wavelength
refl3dat$source <- "Option 3"
#head(refl3dat)

refl4dat <- data.frame(Reflectance=t(refl4dat.m[-1]),sd=t(refl4dat.sd[-1]))
refl4dat$Bandnb <- FX17bands$Bandnb
refl4dat$Wavelength <- FX17bands$Wavelength
refl4dat$source <- "Option 4"
#head(refl4dat)

refldat <- rbind(refl1dat, refl2dat, refl3dat, refl4dat)
#head(refldat)
a <- data.frame(Reflectance=RefMCAA01$Reflectance, sd=0.0, 
                 Bandnb=NA, Wavelength=RefMCAA01$Wavelength, source="Certified")
refldat <- rbind(refldat, a)
rm(a)
#micolor <- c( "orange", "red", "blue", "cyan")
#Default colors
ggplotColours <- function(n = 6, h = c(0, 360) + 15){
  if ((diff(h) %% 360) < 1) h[2] <- h[2] - 360/n
  hcl(h = (seq(h[1], h[2], length = n)), c = 100, l = 65)
}

micolor <- ggplotColours(n=4)
micolor <- c("Option 1"=micolor[1], "Option 2"=micolor[2], 
             "Option 3"=micolor[3], "Option 4"=micolor[4], 
             "Certified"= "black" )
ggplot(data=refldat) +
  #geom_line(data=RefMCAA01, aes(x=Wavelength, y=Reflectance), color="black") +
  #geom_point(data=RefMCAA01, aes(x=Wavelength, y=Reflectance), color="black",size=.25) +
  geom_point(aes(x=Wavelength, y=Reflectance, color=source),size=1) +
  geom_line(aes(x=Wavelength, y=Reflectance, color=source)) +
  geom_errorbar(aes(x=Wavelength,ymin=Reflectance-sd, ymax=Reflectance+sd, color=source),
                width=.1,
                position=position_dodge(0.05)) +
  geom_rect(aes(xmin = 1450, xmax = 1600, ymin = 0.4, ymax =0.9),col="red",alpha=0) +
  xlim(range(FX17bands$Wavelength)) + xlab("Wavelength (nm)")+
  theme(legend.position = c(0.91, 0.23)) +
  scale_color_manual(values=micolor) +
  ggtitle("Reference WCS-MC-020 spectra")

ggplot(data=refldat) +
  #geom_line(data=RefMCAA01, aes(x=Wavelength, y=Reflectance), color="black") +
  #geom_point(data=RefMCAA01, aes(x=Wavelength, y=Reflectance), color="black",size=.5) +
  geom_point(aes(x=Wavelength, y=Reflectance, color=source),size=1) +
  geom_line(aes(x=Wavelength, y=Reflectance, color=source)) +
  geom_errorbar(aes(x=Wavelength,ymin=Reflectance-sd, ymax=Reflectance+sd, color=source),
                width=.1,
                position=position_dodge(0.05)) +
  #geom_rect(aes(xmin = 1450, xmax = 1600, ymin = 0.4, ymax =0.9),col="red",alpha=0) +
  xlim(c(1450,1600)) + ylim(c(0.4, 0.9)) + xlab("Wavelength (nm)") +
  theme(legend.position = c(0.85, 0.22)) +
  scale_color_manual(values=micolor) +
  ggtitle("Reference WCS-MC-020 spectra (Inset)")

Reminder: For the calculation of the Reflectance image

• Option 1. Dark noise not taken into account. Illumination heterogeneity not taken into account.
• Option 2. Dark noise taken into account. Illumination heterogeneity not taken into account.
• Option 3. Dark noise not taken into account. Illumination heterogeneity taken into account.
• Option 4. Dark noise taken into account. Illumination heterogeneity taken into account.

From these graphics, we can conclude:

1. The standard deviation is higher in the first two options, which indicates that the variability introduced by the heterogeneous illumination has an impact even within a relatively small object.
2. As the White image has lower values at the extremes, the global average is low and thus the calculated Reflectance is higher for the first two options.
3. Taking into account both the illumination heterogeneity and the dark noise, Option 4 is the closest to the certified spectrum. Note this is the option we are regularly using to calculate Reflectance from images acquired by FX17 + Lumo scanner in HSILab.
4. Note that the band width (fwhm) of the FX17 instrument implies that narrow absorption features in the certified spectrum appear smoother in FX17 spectra.