## ----include = FALSE----------------------------------------------------------
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>",
  warning = FALSE,
  fig_width = 6,
  fig_height = 4,
  dpi = 72
)

data.table::setDTthreads(2)

## ----setup--------------------------------------------------------------------
library(OpenSpecy)

## ----eval=FALSE---------------------------------------------------------------
# run_app()

## ----eval=FALSE---------------------------------------------------------------
# spectra <- read_any("path/to/your/data")

## ----eval=FALSE---------------------------------------------------------------
# cache_dir <- file.path(tempdir(), "openspecy-large-map-cache")
# map <- open_specs("D:/spectral-data/map.h5", cache_dir = cache_dir)
# map
# 
# # Natively copy one region at a time as map_Region1.h5, etc. Spectral values
# # are not loaded into R.
# region_files <- split_h5(
#   "D:/spectral-data/map.h5",
#   output_dir = "D:/spectral-data/regions"
# )
# 
# # These are descriptor-only regional views.
# region_views <- split_spec(map, by = "region")
# 
# # Materialize only a reviewed selection. A whole-source request is refused.
# pixel <- decompress_spec(region_views[[1]], index = 1)
# roi <- decompress_spec(map, region = "Region1",
#                        roi = c(xmin = 10, xmax = 30,
#                                ymin = 20, ymax = 40))
# 
# # Save the descriptor, not a copy of the spectral cube.
# write_specs(map, "map.filespec.rds")
# map <- read_specs("map.filespec.rds")
# 
# # Or stream one complete rectangular region to a new float64 ENVI pair.
# # Existing output members are refused rather than overwritten.
# write_spec(region_views[[1]], "region-1-export.hdr")

## ----eval=FALSE---------------------------------------------------------------
# particle_result <- automate_particle_analysis(
#   map,
#   library = reference_library,
#   file_processing = "stream",
#   particle_id_strategy = "collapse",
#   collapse_function = mean,
#   spectral_smooth = FALSE,
#   sn_threshold_min = 0.04,
#   sn_threshold_max = Inf,
#   area_threshold = 1,
#   cor_threshold = 0.7,
#   outputs = c("details", "summary", "processed", "particle_image",
#               "particle_heatmap", "particle_heatmap_thresholded",
#               "sn_histogram", "cor_histogram")
# )
# 
# # Recorded plots survive after the pipeline returns.
# plot(particle_result, sample = "Region1", which = "particle_heatmap")
# plot(particle_result, sample = "Region1", which = "sn_histogram")
# plot(particle_result, sample = "Region1", which = "cor_histogram")

## -----------------------------------------------------------------------------
data("raman_hdpe")

## -----------------------------------------------------------------------------
spectral_map <- read_extdata("CA_tiny_map.zip") |> 
  read_any() # preserves some metadata
asp_example <- read_extdata("ftir_ldpe_soil.asp") |> 
  read_any() 
ps_example <- read_extdata("ftir_ps.0") |> 
  read_any() # preserves some metadata
csv_example <- read_extdata("raman_hdpe.csv") |> 
  read_any()
json_example <- read_extdata("raman_hdpe.json") |> 
  read_any() # read in exactly as an OpenSpecy object

## -----------------------------------------------------------------------------
scratch_OpenSpecy <- as_OpenSpecy(x = seq(1000,2000, by = 5), 
                                  spectra =  data.frame(runif(n = 201)), 
                                  metadata = list(file_name = "fake_noise")) 

## -----------------------------------------------------------------------------
# Access the wavenumbers
scratch_OpenSpecy$wavenumber

## -----------------------------------------------------------------------------
# Access the spectra
scratch_OpenSpecy$spectra

## -----------------------------------------------------------------------------
# Access the metadata
scratch_OpenSpecy$metadata

## -----------------------------------------------------------------------------
# Performs checks to ensure that OpenSpecy objects are adhering to our standards;
# returns TRUE if it passes. 
check_OpenSpecy(scratch_OpenSpecy)

# Checks only the object type to make sure it has OpenSpecy type
is_OpenSpecy(scratch_OpenSpecy)

## -----------------------------------------------------------------------------
print(scratch_OpenSpecy) # shows the raw object

## -----------------------------------------------------------------------------
summary(scratch_OpenSpecy) # summarizes the contents of the spectra

## -----------------------------------------------------------------------------
head(scratch_OpenSpecy) # shows the top wavenumbers and intensities

## ----eval=FALSE---------------------------------------------------------------
# write_spec(scratch_OpenSpecy, "test_scratch_OpenSpecy.json", digits = 5)
# write_spec(scratch_OpenSpecy, "test_scratch_OpenSpecy.csv", digits = 5)

## ----eval=FALSE---------------------------------------------------------------
# hyperspecy <- as_hyperSpec(scratch_OpenSpecy)

## ----fig.align="center", fig.width=5------------------------------------------
plot(scratch_OpenSpecy) # quick and efficient

## ----eval=FALSE---------------------------------------------------------------
# # This will min-max normalize your data even if it isn't already but are not
# # changing your underlying data
# plotly_spec(scratch_OpenSpecy, json_example)

## ----eval=FALSE---------------------------------------------------------------
# heatmap_spec(spectral_map,
#              z = spectral_map$metadata$x)

## ----eval=FALSE---------------------------------------------------------------
# interactive_plot(spectral_map, select = 100, z = spectral_map$metadata$x)

## -----------------------------------------------------------------------------
combined <- c_spec(list(asp_example, ps_example), range = "common", res = 8) 

combined |> 
  plot()

## -----------------------------------------------------------------------------
plot(combined |> make_rel(), offset = 3, legend_var = "file_name")

## ----eval=FALSE---------------------------------------------------------------
# # Extract the 150th spectrum by index number.
# filter_spec(spectral_map, 150)
# # Extract the 150th spectrum by spectrum name.
# filter_spec(spectral_map, "9_5")
# # Extract the 150th spectrum by logical test.
# filter_spec(spectral_map, spectral_map$metadata$x == 9 &  spectral_map$metadata$y == 5)
# 
# #Test that they are the same.
# identical(filter_spec(spectral_map, 150), filter_spec(spectral_map, "9_5"))
# 
# identical(filter_spec(spectral_map, 150), filter_spec(spectral_map, spectral_map$metadata$y == 9 & spectral_map$metadata$x == 5))

## -----------------------------------------------------------------------------
sample_spec(spectral_map, size = 5) |>
  plot()

## ----eval=FALSE---------------------------------------------------------------
# processed <- process_spec(raman_hdpe)

## ----eval=FALSE---------------------------------------------------------------
# quality_ready <- raman_hdpe |>
#   correct_spike(direction = "both") |>
#   restrict_range(saturation = "auto", make_rel = FALSE) |>
#   process_spec(correct_spike = FALSE)

## ----eval=FALSE---------------------------------------------------------------
# assessment <- assess_spec(
#   processed,
#   checks = c("co2_region", "high_tail")
# )
# corrected <- processed |>
#   flatten_range(automate = TRUE, make_rel = FALSE) |>
#   restrict_range(automate = TRUE, make_rel = FALSE)

## ----eval=FALSE---------------------------------------------------------------
# quality_report <- assess_spec(
#   quality_ready,
#   checks = c("missing_values", "negative_intensity", "low_snr",
#              "spike", "saturation"),
#   report = "all",
#   snr_metric = "breakpoint_snr"
# )

## ----eval=FALSE---------------------------------------------------------------
# plotly_spec(raman_hdpe, processed)

## ----eval=FALSE---------------------------------------------------------------
# # Automatic signal to noise ratio comparison
# sig_noise(processed, metric = "run_sig_over_noise") >
#   sig_noise(raman_hdpe, metric = "run_sig_over_noise")
# 
# #Manual signal to noise ratio calculation
# sig_noise(processed, metric = "sig_over_noise", sig_min = 2700, sig_max = 3000, noise_min = 1500, noise_max = 2500) >
#   sig_noise(raman_hdpe, metric = "sig_over_noise", sig_min = 2700, sig_max = 3000, noise_min = 1500, noise_max = 2500)

## ----eval=FALSE---------------------------------------------------------------
# #Remove CO2 region
# spectral_map_p <- spectral_map |>
#   process_spec(flatten_range = TRUE)
# 
# #Calculate signal times noise
# spectral_map_p$metadata$sig_noise <- sig_noise(spectral_map_p, metric = "run_sig_over_noise")
# 
# #Plot result
# heatmap_spec(spectral_map_p, sn = spectral_map_p$metadata$sig_noise, min_sn = 5)

## ----eval=FALSE---------------------------------------------------------------
# trans_raman_hdpe <- raman_hdpe
# trans_raman_hdpe$spectra <- 2 - trans_raman_hdpe$spectra^2
# 
# rev_trans_raman_hdpe <- trans_raman_hdpe |>
#   adj_intens(type = "transmittance")
# 
# plotly_spec(trans_raman_hdpe, rev_trans_raman_hdpe)

## ----eval=FALSE---------------------------------------------------------------
# conform_spec(raman_hdpe, res = 8) |> # Convert res to 8 wavenumbers.
#   summary()
# 
# # Force one spectrum to have the exact same wavenumbers as another
# conform_spec(asp_example, range = ps_example$wavenumber, res = NULL) |>
#   summary()
# 

## ----smooth_intens, fig.cap = "Sample `raman_hdpe` spectrum with different smoothing polynomials."----
none <- make_rel(raman_hdpe)
p1 <- smooth_intens(raman_hdpe, polynomial = 1, derivative = 0, abs = FALSE)
p4 <- smooth_intens(raman_hdpe, polynomial = 4, derivative = 0, abs = FALSE)

c_spec(list(none, p1, p4)) |> 
  plot()

## ----compare_derivative, fig.cap = "Sample `raman_hdpe` spectrum with different derivatives."----
none <- make_rel(raman_hdpe)
window <- calc_window_points(raman_hdpe, 100) #Calculate the number of points needed for a 190 wavenumber window.
d1 <- smooth_intens(raman_hdpe, derivative = 1, window = window, abs = TRUE)
d2 <- smooth_intens(raman_hdpe, derivative = 2, window = window, abs = TRUE)

c_spec(list(none, d1, d2)) |> 
  plot()

## ----subtr_baseline, fig.cap = "Sample `raman_hdpe` spectrum with different degrees of background subtraction (Cowger et al., 2020)."----
alternative_baseline <- smooth_intens(raman_hdpe, polynomial = 1, window = 51,
                                      derivative = 0, abs = FALSE,
                                      make_rel = FALSE) |>
  flatten_range(min = 2700, max = 3200, make_rel = FALSE) #Manual baseline with heavily smoothed spectra

none <- make_rel(raman_hdpe) #raw
d <- subtr_baseline(raman_hdpe, type = "manual",
                     baseline = alternative_baseline) #manual subtraction
d8 <- subtr_baseline(raman_hdpe, degree = 8) #standard imodpolyfit
dr <- subtr_baseline(raman_hdpe, refit_at_end = TRUE) #optionally retain baseline noise with refitting

c_spec(list(none, d, d8, dr)) |> 
  plot(offset = 0.25)

## -----------------------------------------------------------------------------
fill_peaks <- subtr_baseline(
  asp_example,
  type = "fill_peaks", lambda = 4, hwi = 50, it = 10,
  make_rel = FALSE
)

## ----restrict_range, fig.cap = "Sample `raman_hdpe` spectrum with different degrees of range restriction."----
none <- make_rel(raman_hdpe)

#Specify one range
r1 <- restrict_range(raman_hdpe, min = 1000, max = 2000) |>
    conform_spec(range = none$wavenumber, res = NULL, allow_na = TRUE)

#Specify multiple ranges
r2 <- restrict_range(raman_hdpe, min = c(1000, 1800), max = c(1200, 2000)) |>
        conform_spec(range = none$wavenumber, res = NULL, allow_na = TRUE)

compare_ranges <- c_spec(list(none, r1, r2), range = "common")
# Common argument crops the ranges to the most common range between the spectra
# when joining. 

plot(compare_ranges)

## ----flatten_range, fig.cap = "Sample `raman_hdpe` spectrum with different degrees of background subtraction (Cowger et al., 2020)."----
single <- filter_spec(spectral_map, 120) # Function to filter spectra by index
                                         # number or name or a logical vector. 
none <- make_rel(single)
f1 <- flatten_range(single) #default flattening the CO2 region. 
f2 <- flatten_range(single, min = c(1000, 2500), max = c(1200, 3000)) #multple range example

compare_flats <- c_spec(list(none, f1, f2), range = "common")

plot(compare_flats, offset = 0.25)

## ----eval=FALSE---------------------------------------------------------------
# 
# raman_hdpe |> plot()
# 
# make_rel(raman_hdpe) |> plot()
# 

## -----------------------------------------------------------------------------
single_region_area <- area_under_band(
  fill_peaks, min = 2700, max = 3000
)
single_point_intensity <- point_intensity(
  fill_peaks, wavenumber = 2880
)
single_region_area
single_point_intensity

## -----------------------------------------------------------------------------
custom_area_numerator <- area_under_band(
  fill_peaks, min = 2700, max = 3000
)
custom_area_denominator <- area_under_band(
  fill_peaks, min = 1400, max = 1500
)
custom_area_ratio <- custom_area_numerator / custom_area_denominator
custom_area_ratio[!is.finite(custom_area_ratio)] <- NA_real_
custom_area_ratio

## -----------------------------------------------------------------------------
saub_ci <- area_under_band(
  fill_peaks,
  index = "carbonyl_index_saub"
)
saub_ci

## -----------------------------------------------------------------------------
nearest_point_ratio <- peak_ratio(
  fill_peaks, numerator = 2880, denominator = 2840
)
linear_point_ratio <- peak_ratio(
  fill_peaks, numerator = 2880, denominator = 2840, method = "linear"
)
nearest_point_ratio
linear_point_ratio

## ----eval=FALSE---------------------------------------------------------------
# get_lib(type = "derivative")

## ----eval=FALSE---------------------------------------------------------------
# lib <- load_lib(type = "derivative")

## ----eval = FALSE-------------------------------------------------------------
# data("test_lib")
# data("raman_hdpe")
# 
# processed <- process_spec(x = raman_hdpe,
#                           conform_spec = FALSE, #We will conform during matching.
#                           smooth_intens = TRUE #Conducts the default derivative transformation.
#                           )
# 
# # Check to make sure that the signal to noise ratio of the processed spectra is
# # greater than 10.
# print(sig_noise(processed) > 10)
# 
# #Plot to assess the accuracy of the processing visually
# plotly_spec(raman_hdpe, processed)

## ----eval=FALSE---------------------------------------------------------------
# matches <- match_spec(x = processed, library = test_lib, conform = TRUE,
#                       add_library_metadata = "sample_name", top_n = 1,
#                       top_n_by = "organization")[order(match_val, decreasing = TRUE)]
# print(matches[,c("object_id", "library_id", "match_val", "SpectrumType",
#                  "SpectrumIdentity")])

## ----eval=FALSE---------------------------------------------------------------
# get_metadata(x = test_lib, logic = matches[[1,"library_id"]], rm_empty = TRUE)

## ----eval=FALSE---------------------------------------------------------------
# plotly_spec(processed, filter_spec(test_lib, logic = matches[[1,"library_id"]]))

## ----eval = FALSE-------------------------------------------------------------
# #Test library
# data("test_lib")
# 
# #Example hyperspectral image with one cellulose acetate particle in the middle of it.
# test_map <- read_any(read_extdata("CA_tiny_map.zip"))
# 
# #Process the map to conform to the library.
# test_map_processed <- process_spec(test_map, conform_spec_args = list(
#   range = test_lib$wavenumber, res = NULL)
#   )
# 
# #Identify every spectrum in the map.
# identities <- match_spec(test_map_processed, test_lib, order = test_map,
#                          add_library_metadata = "sample_name", top_n = 1)
# 
# #Relabel any spectra with low correlation coefficients.
# features <- ifelse(identities$match_val > 0.7,
#                    tolower(identities$polymer_class), "unknown")
# 
# #Use spectra identities to identify particle regions as those that have the same material type and are touching.
# id_map <- def_features(x = test_map_processed, features = features)
# 
# id_map$metadata$identities <- features
# # Also should probably be implemented automatically in the function when a
# # character value is provided.
# heatmap_spec(id_map, z = id_map$metadata$identities)
# 
# # Collapses spectra to their median for each particle
# test_collapsed <- collapse_spec(id_map)
# 
# # Plot spectra for each identified particle
# plot(test_collapsed, offset = 1, legend_var = "feature_id")

## ----eval = FALSE-------------------------------------------------------------
# # Read in test library
# data("test_lib")
# 
# # Example dataset with one cellulose acetate particle.Conduct spatial smoothing to average each spectrum using adjacent spectra.
# test_map <- read_any(read_extdata("CA_tiny_map.zip"),
#                      spectral_smooth = TRUE,
#                      sigma = c(1, 1, 1))
# 
# # Characterize the signal times noise to determine where particle regions are.
# snr <- sig_noise(test_map, metric = "sig_times_noise")
# 
# # Use this to find your particles and the idal signal times noise value to use for thresholding.
# heatmap_spec(test_map, z = snr)
# 
# # Define the feature regions based on the threshold. Pixels from the background in the heatmap above were below 0.05 while my particle's pixels were above so I set snr > 0.05.
# id_map <- def_features(x = test_map, features = snr > 0.05)
# 
# # Check that the thresholding worked as expected. Here we see a single particle region identified separate from the background.
# heatmap_spec(id_map, z = id_map$metadata$feature_id)
# 
# # Collapse the spectra to their medians based on the threshold. Important to
# # note here that the particles with id -88 are anything from the FALSE values
# # so they should be your background.
# collapsed_id_map <- id_map |>
#   collapse_spec()
# 
# # Process the collapsed spectra to have the same transformation and units as the library.
# id_map_processed <- process_spec(collapsed_id_map, conform_spec_args = list(
#   range = test_lib$wavenumber, res = NULL)
#   )
# 
# # Check the spectra for the background and particle. Background has considerable signal in it too suggesting double bounce along the edges of the particle.
# plot(id_map_processed, offset = 1, legend_var = "feature_id")
# 
# # Get the matches of the collapsed spectra for the particles.
# matches <- match_spec(id_map_processed, test_lib,
#                       add_library_metadata = "sample_name", top_n = 1)
#