---
title: "Open Specy Package Tutorial"
author: >
  Win Cowger, Zacharias Steinmetz, Rachel Kozloski, Aleksandra Karapetrova
date: "`r Sys.Date()`"
output: rmarkdown::html_vignette
vignette: >
  %\VignetteIndexEntry{Open Specy Package Tutorial}
  %\VignetteEncoding{UTF-8}
  %\VignetteEngine{knitr::rmarkdown}
---

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

data.table::setDTthreads(2)
```

# Document Overview 

This document outlines a common workflow for using the Open Specy package and
highlights some topics that users are often requesting a tutorial on. If the
document is followed sequentially from beginning to end, the user will have a
better understanding of every procedure involved in using the Open Specy R
package as a tool for interpreting spectra. It takes approximately 45 minutes to
read through and follow along with this standard operating procedure the first
time. Afterward, knowledgeable users should be able to thoroughly analyze
spectra at an average speed of 1 min^-1^ or faster with the new batch and
automated procedures.

The Open Specy R package is the backbone of the Shiny app. The choice is yours
as to which you start with, we use both on a regular basis. The tutorial will
talk through the R functions you can use to programatically analyze spectra. If
you are looking for a tutorial about how to use the app see [this
tutorial](https://www.openanalysis.org/OpenSpecyV2/pkgdown/articles/app.html).

# Installation

After installing the R package, you just need to read in the library.

```{r setup}
library(OpenSpecy)
```

# Running the App

To get started with the Open Specy user interface, access
[the embedded browser app](https://www.openanalysis.org/OpenSpecyV2/) or
start the Shiny GUI directly from your own computer in R. If you are looking for
a tutorial about how to use the app see [this
tutorial](https://www.openanalysis.org/OpenSpecyV2/pkgdown/articles/app.html).

```{r, eval=FALSE}
run_app()
```

# Read Data

The following line of code will read in your data when using the package and
interprets which reading function to use based on the file extension.

```{r, eval=FALSE}
spectra <- read_any("path/to/your/data")
```

Open Specy allows for upload of native Open Specy .csv, .json, or .rds files. Open Specy and .csv files should always load correctly
but the other file types are still in development, though most of the time these
files work perfectly.
In addition, .csv, .asp, .jdx, .0, .spa, .spc, and .zip files can be imported.
.zip  files can either contain multiple files with individual spectra in them of
the non-zip formats or it can contain a .hdr and .dat file that form an ENVI
file for a spectral map. If uploading a .csv file, it is ideal to label the column with the
wavenumbers `wavenumber` and name the column with the intensities `intensity`. Columns besides wavenumber will be interpreted as unique spectra. If any columns are numbers, the csv will be interpreted in wide format with the number columns the wavenumbers and rows containing the unique spectral intensities and metadata contained in non number columns.
Wavenumber units should be cm^-1^ or use `adj_wave()` to correct from wavelength. 
Always keep a copy of the original file before alteration to preserve metadata
and raw data for your records.

It is best practice to cross check files in the proprietary software they came
from and Open Specy before use in Open Specy. Due to the complexity of some
proprietary file types, we haven't been able to make them fully compatible yet.
If your file is not working, please contact the administrator and share the file
so that we can work on integrating it.

The specific steps to converting your instrument's native files to .csv can be
found in its software manual or you can check out
[Spectragryph](https://www.effemm2.de/spectragryph/), which supports many
spectral file conversions. For instructions, see [Spectragryph
Tutorial](https://www.openanalysis.org/OpenSpecyV2/pkgdown/articles/spectragryph.html). Unfortunately the maintainer of Spectragryph passed away and it is unclear how much longer this will be supported.

## Experimental large H5 and ENVI workflow

Use `FileSpecs` for a hyperspectral map that is too large for an ordinary
in-memory `OpenSpecy` workflow. It is a lightweight, local-first `Specs`
subtype: the serialized object contains source, axis, coordinate, recipe,
fingerprint, and cache descriptors, but no open connection and no full values
matrix. `open_specs()` currently accepts H5 or an ENVI `.hdr` plus `.dat`/`.img`
pair. It opens each member read-only and never writes beside or over the source.

Choose a dedicated cache directory. Completed cache entries are versioned by
the source fingerprint and published atomically. A changed, missing, or moved
source is rejected instead of silently reusing stale results. A saved
descriptor therefore still needs its original, unchanged source files. Remove
only a cache directory that you created, and only after every session using it
has ended; source files are not cache files.

```{r, 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")
```

`split_h5()` defaults to `field = "region"` and native H5 output. Another
indexed or H5 metadata field can be used when each category contains complete
regions. For a field that divides pixels within a region, use `format = "rds"`;
each materialized category must then fit in memory. The function returns output
paths invisibly, sanitizes invalid filename characters, and refuses to
overwrite existing files. Native splits retain the registered H5 mosaic tiles
that intersect their selected regions, with matching `Centers` rows and
renumbered image datasets, so `visual_image()` and particle overlays continue
to work on each output.

The export preserves the wavelength axis with round-trip-safe formatting and
rejects custom writers, incomplete/multi-region views, source members, and
existing output paths.

Selecting an entire region is explicit but may still be large; prefer an index
or a small region of interest before calling established `OpenSpecy` functions.
Matrix-only `Specs` operations fail with guidance rather than falling back to a
hidden whole-cube read.

The first supported direct analysis is region-aware particle identification
with S/N screening and exact mean collapse. It streams each region, rereads
only retained pixels to form particle means, materializes that small collapsed
object, and then uses established preprocessing and exact matching. Thresholds
remain scientific choices: `sn_threshold_min` and `sn_threshold_max` bound the
accepted S/N scale, `area_threshold` rejects smaller connected regions in
pixels, and `cor_threshold` is the minimum correlation used for interpretation.

```{r, 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")
```

For a `FileSpecs` object or H5/HDF5 path, `file_processing = "stream"` is the
default. It reports region, chunk, spectrum-count, percentage, and elapsed-time
milestones while reading bounded blocks. Choose `file_processing = "memory"`
to materialize every region before analysis when the full spectra and analysis
intermediates fit comfortably in RAM. Memory mode enables the complete
in-memory strategy set and can be faster for smaller maps, but it removes the
streaming memory bound.

The signal heatmap uses a continuous finite minimum-to-maximum legend; the
threshold mask remains binary. Both histograms draw the configured threshold
lines. The pipeline retains the exact best match on each collapsed particle.
For H5 data, requesting
`particle_image` materializes and caches the registered visual mosaic only for
the current region, using all intersecting tiles rather than one unrelated
tile. Stream mode supports connected `collapse` with `mean` and
`all_cell_id`, non-entropy S/N, and bounded spectral smoothing. Other
strategies, median/custom collapse, and arbitrary whole-map operation chains
require memory mode. Verify stream/memory equivalence on a representative
subset before interpreting a large run.

For an ENVI DAT/IMG path, `automate_particle_analysis()` also looks beside the
spectral pair for one same-basename JPG, JPEG, or PNG when `images` is omitted.
It uses the red-boundary detector to register that image and draws the same
material-class overlay as `particle_image()`. Multiple matches or an
undetectable boundary warn and leave the spectral analysis usable without an
image.

Local Shiny uses one **Choose spectra...** control: a native Windows/macOS
multi-file dialog when available, or a mounted-filesystem browser elsewhere.
It keeps the original paths instead of
uploading duplicate files. Shinylive mounts browser-selected files through
WORKERFS. Its 10 GiB input ceiling is a transport limit rather than a promise
that every operation will fit in RAM. ENVI/H5 maps normally use compact
file-backed coordinates and metadata; **Load Entire File into Memory** is an
explicit advanced override, with a configurable identification batch size.
Optional background suppression maps rejected pixels to virtual zero while
retaining per-source S/N and reasons. Selection only stages a source; **Run**
owns the read and every later operation.

For an uploaded map, the app applies optional spatial smoothing before any
particle decision. Signal/noise is calculated from that spatial-only object.
Connected regions collapse eligible pixels before ordinary processing and
identification. Both clustering strategies instead fit one source-scoped PCA
plus K-means model to the spatial-only spectra, collapse those spectral groups,
and process and identify the clusters once. Non-spatial mode returns those
clusters. Spatial mode projects the cluster material identities to the original
pixels, joins adjacent equal-material clusters, collapses the raw/spatial-only
spectra again, and reprocesses the final particles without identifying twice.
When correlation is a threshold, that same first Top N identification result is
reused. Threshold failures are background pixels, not unknown particles.

The default **Mean Up** conformation technique leaves the uploaded axis
unchanged and conforms the identification library onto it with memory-bounded
averaging and interpolation, unless the selected Wavenumber Resolution is
finer than the upload's native resolution, in which case the uploaded spectra
are resampled up to that resolution instead. Heatmap legends open from **View Legend**; categorical
maps with more than 30 values rely on hover instead of listing an unreadable
legend. A rejected pixel has no nearest-particle substitution: its processed
trace is flat and it has no matches. Thresholded-particle downloads can include
the final summary table, both histograms, all available heatmaps, the material
summary, and particle-size distribution. The package-only experimental workflow
above remains available independently for bounded on-disk research.

Maintainers can run the repeated benchmark in
`benchmarks/file_specs_particle_analysis.R`. It is opt-in and accepts external
H5, ENVI, eager-RDS oracle, and reference-library paths through environment
variables; it never copies the large source into the repository.

If you don't have your own data, you can use a test dataset.

```{r}
data("raman_hdpe")
```

We also have many onboard files that you can call to test different formats:

```{r}
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
```

You will notice now that the R package reads in files into an object with class
`OpenSpecy`. This is a class we created for high throughput spectral analysis
which now also preserves spectral metadata. You can even create these from
scratch if you'd like.

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

Open Specy objects are lists with three components, `wavenumber` is a vector of
the wavenumber values for the spectra and corresponds to the rows in `spectra`
which is a `data.table` where each column is a set of spectral intensities.
`metadata` is a data.table which holds additional information about the spectra.
Each row in `metadata` corresponds to a column in `spectra`.

```{r}
# Access the wavenumbers
scratch_OpenSpecy$wavenumber
```

```{r}
# Access the spectra
scratch_OpenSpecy$spectra
```

```{r}
# Access the metadata
scratch_OpenSpecy$metadata
```

```{r}
# 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)
```

We have some generic functions built for inspecting the spectra:

```{r}
print(scratch_OpenSpecy) # shows the raw object
```

```{r}
summary(scratch_OpenSpecy) # summarizes the contents of the spectra
```

```{r}
head(scratch_OpenSpecy) # shows the top wavenumbers and intensities
```

# Save Data

Open Specy objects can be saved most accurately as .csv, .rds, or .json files.
.rds will be the most reproducible as it is a native R file format and floating
point errors can happen with .csv or .json.

```{r, eval=FALSE}
write_spec(scratch_OpenSpecy, "test_scratch_OpenSpecy.json", digits = 5)
write_spec(scratch_OpenSpecy, "test_scratch_OpenSpecy.csv", digits = 5)
```

# Format Conversions

Another great spectroscopy R package is
[hyperSpec](https://CRAN.R-project.org/package=hyperSpec). We
actually depend on their functions for several of ours. They are currently
making some awesome new features and we want to integrate well with them so that
both packages can be easily used together. That is why we created the
`as_hyperSpec` function and made sure that `as_OpenSpecy` can convert from
`hyperSpec` objects.

```{r, eval=FALSE}
hyperspecy <- as_hyperSpec(scratch_OpenSpecy)
```

# Visualization

## Spectra

In R, we have two ways to visualize your spectra, one is quick and efficient and
the other is interactive. Here is an example of quick and efficient plotting.

```{r, fig.align="center", fig.width=5}
plot(scratch_OpenSpecy) # quick and efficient
```

This is an example of an interactive plot. You can plot two different datasets
simultaneously to compare. 

```{r, 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)
```

## Maps

Spectral maps can also be visualized as overlaid spectra but in addition the
spatial information can be plotted as a heatmap. It is important to
note that when multiple spectra are uploaded in batch they are prescribed `x`
and `y` coordinates, this can be helpful for visualizing summary statistics and
navigating vast amounts of data but shouldn't be confused with data which
actually has spatial coordinates.

In R the user can control what values form the colors of the heatmap by
specifying `z`, it is handy to put the information you want in the metadata for
this reason. This example just shows the x values of the spectra.

```{r, eval=FALSE}
heatmap_spec(spectral_map,
             z = spectral_map$metadata$x)
```

An interactive plot of the heatmap and spectra overlayed can be generated with
the `interactive_plot()` function. A user can hover over the heatmap to identify
the next row id they are interested in and update the value of select to see
that spectrum.

```{r, eval=FALSE}
interactive_plot(spectral_map, select = 100, z = spectral_map$metadata$x)
```
<br>

# Combining OpenSpecy Objects

Sometimes you have several OpenSpecy objects that you want to combine into one.
The easiest way to do that is by having spectra which all are in the same exact
format with the same series of wavenumbers. The default settings of `c_spec`
assume that is the case.

If you have different wavenumber ranges for the spectra you want to combine,
you can set  `range = "common"` and `res` equal to the wavenumber resolution you
want and the function will collapse all the spectra to whatever their common
range is using linear interpolation.

```{r}
combined <- c_spec(list(asp_example, ps_example), range = "common", res = 8) 

combined |> 
  plot()
```

If you want to plot them with standardized offset y values you can do so with the plot variables. 

```{r}
plot(combined |> make_rel(), offset = 3, legend_var = "file_name")
```


# Filtering OpenSpecy Objects

OpenSpecy objects can have any number of spectra in them. To create an OpenSpecy
with a subset of the spectra that is in an Open Specy object you can use the
`filter_spec` function. Filtering is allowed by index number, name, or using a
logical vector. Filtering will update the `spectra` and `metadata` items of the
`OpenSpecy` but not the `wavenumber`.

```{r, 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))
```

# Sampling OpenSpecy Objects

Sometimes you have a large suite of examples of spectra of the same material and
you want to reduce the number of spectra you run through the analysis for
simplicity or you are running simulations or other procedures that require you
to first sample from the spectra contained in your OpenSpecy objects before
doing analysis. The `sample_spec` function serves this purpose.

```{r}
sample_spec(spectral_map, size = 5) |>
  plot()
```

# Processing

The goal of processing is to increase the signal to noise ratio (S/N)
of the spectra and remove unwanted artifacts without distorting the shape, position, or relative size of the
peaks. After loading data, you can process the data using intensity
adjustment, baseline subtraction, smoothing, flattening, and range selection.
The default settings is an absolute first derivative transformation. It is really powerful for many data issues. It does something similar to smoothing, baseline subtraction, and
intensity correction simultaneously and really quickly.

The `process_spec()` function is a monolithic function for all
processing procedures which is optimized by default to result in high signal to
noise in most cases, same as the app.

```{r eval=FALSE}
processed <- process_spec(raman_hdpe)
```

Isolated impulse artifacts and detector saturation should be handled before
ordinary smoothing, baseline correction, derivatives, or normalization. The
package keeps both operations opt-in; the app exposes them as owner controls and
uses this order when they are enabled:

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

The default method, contributed by Nicolas Coca Lopez, detects positive and
negative narrow peaks using prominence against the raw median absolute
deviation of first differences. Its maximum spike width defaults to 2 points,
its noise multiplier to 10, and its interpolation window to 5 points per side.
It is implemented in base R, without `pracma`. The previous conservative,
wavenumber-aware detector remains available as `method = "residual"`.
`correct_spike()` also retains the manual and ratio-based prominence/FWHM
procedures described by
[Coca-Lopez (2024)](https://doi.org/10.1016/j.aca.2024.342312). Manual thresholds
must be calibrated for the material, instrument, resolution, and acquisition
conditions. Width and interpolation windows are measured in points, not cm^-1^.

No single-spectrum detector can always distinguish a cosmic-ray spike from a
genuine narrow band. Inspect correction diagnostics, calibrate on representative
standards, and pay particular attention to spikes on real bands and to
narrow-band materials such as calcite and polystyrene. Boundary candidates are
not extrapolated. If a feature is ambiguous, retain the raw data and recollect
rather than treating an interpolated peak as measured evidence.

`saturation = "auto"` detects a hard maximum plateau only when it has lower
shoulders; a known finite detector ceiling can be supplied instead. Detected
intervals are unioned across the batch, guard-expanded, and removed from every
spectrum on the same shared axis. At most 70% of wavenumber-cell coverage is
removed by default. A proposal over that inclusive limit, or one leaving too
few points, is rejected with a warning and unchanged spectral values. Apply the
identical retained axis to the reference library before identification. Open
Specy does not reconstruct saturated peaks: recollection at lower exposure,
laser power, or integration time is the scientifically preferred remedy.

Processed spectra can be checked and corrected before identification. The
default CO2 and high-tail checks compare the artifact-region maximum with the
maximum outside both artifact regions after temporary min-max normalization.
The default artifact ratio is 2: a candidate is flagged at or above twice the
control-region maximum. Raising this dimensionless ratio makes detection and
automatic correction less sensitive; an explicit value remains authoritative.
The default silent-region quality check covers 2420--2550 cm^-1^ and is applied
to each assessed spectrum, not as a blanket dataset result.
With `automate = TRUE`, each range function returns spectra unchanged when its
issue is absent. Otherwise, `flatten_range()` flattens a flagged CO2 region and
`restrict_range()` crops a shared batch tail only when no more than 20% of the
full wavenumber span must be removed.

```{r 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)
```

The compatibility default returns issue rows only and stays empty when a clean
spectrum has no findings. Use the complete report for auditable UI or batch
status, and opt into spike and saturation checks explicitly. The exact
breakpoint SNR sorts amplitudes and evaluates every possible two-group split;
it does not depend on histogram bins.

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

You can compare the processed and unprocessed data in an overlay plot.

```{r, eval=FALSE}
plotly_spec(raman_hdpe, processed)
```

We want people to use the `process_spec()` function for most processing
operations. All other processing functions can be tuned using its parameters in
the single function see `?process_spec()` for details. However, we recognize that
nesting of functions and order of operations can be useful for users to control
so you can also use individual functions for each operation if you'd like. See
explanations of each processing sub-function below.

## Threshold Signal and Noise

Considering whether you have enough signal to analyze spectra in the first place is important, because if you don't have enough signal you should recollect the spectrum.
Classical spectroscopy would recommend your highest peak to be at least 10 times greater than
the baseline of your processed spectra before you begin analysis. Setting a threshold can assist in standardized differentiation between low and high quality spectral. If your
spectra is below that threshold even after processing, you may want to consider
recollecting it. In practice, we are rarely able to collect spectra of that good
quality and more often use 4 as a lower bound with anything below 2 being completely unusable. The "run_sig_over_noise" `metric` searches your
spectra for high and low regions and conducts division on them to derive the
signal to noise ratio. In the
example below you can see that our signal to noise ratio is increased by the
processing, the goal of processing is generally to maximize the signal to noise ratio. If you know where your signal region and noise regions are, you
can specify them with `sig_min`, `sig_max`, `noise_min`, and `noise_max`.

```{r, 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)
```

If analyzing spectra in batch, we recommend looking at the heatmap and
optimizing the percent of spectra that are above your signal to noise threshold
to determine the correct settings instead of looking through spectra
individually. Setting the `min_sn` will threshold the heatmap image to only
color spectra which have a `sn` value over the threshold.

```{r, 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)
```

## Intensity Adjustment

Most functions in Open Specy assume that intensity units are in absorbance units and Open Specy
can adjust reflectance or transmittance spectra to absorbance units. The
transmittance adjustment uses the $\log_{10} 1/T$ calculation which does not
correct for system or particle characteristics. The reflectance adjustment uses
the Kubelka-Munk equation $\frac{(1-R)^2}{2R}$.

This is the respective R code for a scenario where the spectra doesn't need
intensity adjustment:

```{r, 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)
```

## Conforming

Conforming spectra is essential before comparing to a reference library and can
be useful for summarizing data when you don't need it to be highly resolved
spectrally. We set the default spectral resolution to 5 because this tends to be
pretty good for a lot of applications and is in between 4 and 8 which are
commonly used wavenumber resolutions. 

```{r, 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()

```

## Smoothing

The [Savitzky-Golay
filter](https://en.wikipedia.org/wiki/Savitzky%E2%80%93Golay_filter) is used for
smoothing. Higher polynomial numbers lead to more wiggly fits and thus less
smoothing, lower numbers lead to more smooth fits. The SG filter is fit to a
moving window of 11 data points by default where the center point in the window
is replaced with the polynomial estimate. Larger windows will produce smoother
fits. The derivative order is set to 1 by default which transforms the spectra
to their first derivative. A zero order derivative will have no derivative
transformation and only apply smoothing. When smoothing is done well, peak shapes and relative heights
should not change. The absolute value is primarily useful for first derivative
spectra where the result looks similar to absorbance units and is easy to have intuition about. 

Examples of smoothing:

```{r 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()
```

Derivative transformation can be done with the same function. 

```{r 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()
```

## Baseline Correction

The goal of baseline correction is to get all non-peak regions of the spectra to
zero absorbance. Open Specy provides modified polynomial fitting and 4S Fill
Peaks (Liland 2015) as automated methods. Fill Peaks uses iterative local-window
suppression and is useful for nonlinear baselines with local variation. The
modified polynomial method iteratively fits polynomial equations to the whole
spectrum; a higher polynomial order produces a more flexible fit but can remove
peaks when the baseline is not very wiggly. Manual baseline correction can also
be specified by providing a baseline `OpenSpecy` object. There are many fine
tuning options that can be chosen, see `?subtr_baseline()` for more details.

```{r 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)
```

For example, apply Fill Peaks to the FTIR LDPE spectrum loaded earlier while
retaining its intensity scale for later quantification.

```{r}
fill_peaks <- subtr_baseline(
  asp_example,
  type = "fill_peaks", lambda = 4, hwi = 50, it = 10,
  make_rel = FALSE
)
```

## Range Selection

Sometimes an instrument operates with high noise at the ends of the spectrum
and, a baseline fit produces distortions, or there are regions of interest for
analysis. Range selection accomplishes those goals. Many of these issues can be resolved during spectral collection by specifying a wavenumber range which is well characterized by the instrument. Multiple ranges can
be specified simultaneously.

```{r 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)
```

## Flattening Ranges

Sometimes there are peaks that really shouldn't be in your spectra and can
distort your interpretation of the spectra but you don't necessarily want to
remove the regions from the analysis because you believe those regions should exist and be
flat instead of having a peak. One way to deal with this is to replace the peak
values with the mean of the values around the peak. This is the purpose of the
`flatten_range` function. By default it is set to flatten the CO2 region for
FTIR spectra because that region often needs to be flattened when atmospheric
artifacts occur in spectra. Like `restrict_range`, the R function can accept
multiple ranges.

```{r 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)
```

## Min-Max Normalization

Often we regard spectral intensities as arbitrary and min-max normalization
allows us to view spectra on the same scale without drastically distorting their
shapes or relative peak intensities. In the package, most of the processing
functions will min-max transform your spectra by default if you do not specify
otherwise with `make_rel = FALSE`. 

```{r, eval=FALSE}

raman_hdpe |> plot()

make_rel(raman_hdpe) |> plot()

```

## Area and Point Measurements and Ratios

For a non-ratio measurement, `area_under_band()` returns the inclusive summed
intensity within one explicit region, while `point_intensity()` returns the
intensity at one requested wavenumber for every spectrum. The point method uses
the nearest measured value by default; `method = "linear"` interpolates between
adjacent measured points without extrapolating.

```{r}
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
```

For a custom area-over-area ratio, choose the numerator and denominator ranges
explicitly, calculate each with `area_under_band()`, and divide the returned
vectors. Bounds are inclusive, and each returned value stays aligned with the
corresponding spectrum column. The example below uses arbitrary teaching ranges;
select scientifically appropriate bands and preprocessing for your application.

```{r}
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
```

Named area-under-band ratios select sourced numerator and denominator bands;
they do not apply baseline correction, normalization, or another preprocessing
step. Apply the desired treatment first, then calculate the index. For example,
the specified area-under-band carbonyl index for PE and PP can be calculated
from the Fill Peaks result created above (Almond et al. 2020).

```{r}
saub_ci <- area_under_band(
  fill_peaks,
  index = "carbonyl_index_saub"
)
saub_ci
```

Additional area-over-area presets are `carbonyl_index_pe`,
`hydroxyl_index_pe`, `hydroxyl_index_pp`, `carbon_oxygen_index_pe`, and
`carbon_oxygen_index_pp`. See `?area_under_band` for their sourced band
definitions and intended materials. Baseline treatment and other preprocessing
can materially change an index value, so compare spectra only when the same
treatment was used and keep that treatment with the resulting metadata.

For point ratios, `peak_ratio()` divides the intensity at a numerator
wavenumber by the intensity at a denominator wavenumber for every spectrum.
The default `method = "nearest"` uses measured points and resolves an exact
midpoint tie toward the lower wavenumber. Set `method = "linear"` to interpolate
only between the adjacent measured points. Neither method extrapolates beyond
the available wavenumber axis, and the function does not search for a local
maximum.

```{r}
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
```

# Identifying Spectra 

## Reading Libraries

Reference libraries are spectra with known identities. The Open Specy library
now has more than 59,000 spectra and is too large to bundle in the R package. We
distribute the runtime libraries through AWS and provide `get_lib()` to download
them automatically. Running `get_lib()` by itself downloads every supported
runtime library to your package directory, or you can specify which libraries
you want and where you want them.

```{r, eval=FALSE}
get_lib(type = "derivative")
```

After download you can load the libraries into your active environment one at a time as `OpenSpecy` objects. You
can use any Open Specy object as a library which makes it easy to work with and create libraries because everything we explained earlier applies to them.

```{r, eval=FALSE}
lib <- load_lib(type = "derivative")
```

## Matches

Before attempting to use a reference library to identify spectra it is really
important to understand what format the reference library is in. All the
OpenSpecy reference libraries are in Absorbance units. `derivative` has been
absolute first derivative transformed, `nobaseline` has been baseline corrected,
`raw` is the rawest form of the reference spectra (not recommended except for
advanced uses). The previously mentioned libraries all have Raman and FTIR
spectra in them. `mediod_{derivative or nobasline}` is the mediod compressed library version of the libraries and has only critical spectra in it, `model{derivative or nobasline}` is an exception because it is a
multinomial regression approach for identification. 
In this example we use the `data("test_lib")` which is a subsampled version of
the `derivative` library and `data("raman_hdpe")` which is an unprocessed
Raman spectrum in absorbance units of HDPE plastic. 

```{r, 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)
```

After your spectra is processed similarly to the library specifications, you can
identify the spectra using `match_spec()`. Whichever library you choose, you need to get your
spectra into a similar enough format to use for comparison. The
`add_library_metadata` and `add_object_metadata` options specify the column name
in the metadata that you want to add metadata from and `top_n` specifies how
many matches you want. Set `top_n_by` to one nonblank library metadata column,
such as `"organization"`, to retain that many matches independently from every
group. In this example we just identified a single spectrum with
the library but you can also send an OpenSpecy object with multiple spectra. The
output `matches` is a data.table with at least 3 columns, `object_id` tells you
the column names of the spectra in `x`, `library_id` tells you the column names
from the library that it matched to. `match_val` is the value of the Pearson
correlation coefficient (default) or other correlation if specified in `...` or
if using the model identification option `match_val` will be the model
confidence. The output in this example returned the correct material type, HDPE,
as the top match. If using Pearson correlation, 0.7 is a good threshold to use
for a positive ID. In this example, only our top match is greater than the
threshold so we would disregard the other matches. If no matches were above our
threshold, we would proclaim that the spectrum is of an unknown identity. You'll
also notice in this example that we matched to a library with both Raman and
FTIR spectra but the Raman spectra had the highest hits, this is the rationale
for lazily matching to a library with both. If you want to just match to a
library with FTIR or Raman spectra, you can first filter the library using
`filter_spec()` using `SpectrumType`.

```{r, 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")])
```

## Library Metadata

The libraries we have created have over 100 variables of metadata in them and
this can be onerous to read through especially given that many of the variables
are `NA` values. We created `get_metadata()` to remedy this by removing columns
from the metadata which are all blank values. The function below will return the
metadata for the top match in `matches`. Remember, similar `filter_spec()`, you
can specify `logic` for more than one thing at a time.

```{r, eval=FALSE}
get_metadata(x = test_lib, logic = matches[[1,"library_id"]], rm_empty = TRUE)
```

## Plot Matches

Overlaying unknown spectra and the best matches can be extremely useful to
identify peaks that don't fit to the reference library which may need further
investigation. The example below shows great correspondence between the best
match and the unknown spectrum. All major peaks are accounted for and the
correct relative height. There are two small peaks in the unknown spectrum near
500 that are not accounted for which could be investigated further but we would
call this a positive id to HDPE.

```{r, eval=FALSE}
plotly_spec(processed, filter_spec(test_lib, logic = matches[[1,"library_id"]]))
```

## Sharing Reference Data

If you have reference data or AI models that you think would be useful for
sharing with the spectroscopy community through OpenSpecy please contact the
package administrator to discuss options for collaborating.

# Characterizing Particles

Sometimes the spectroscopy task we want to perform is to identify particles in a
spectral map. This is especially common for microplastic analysis where a
spectral map is used to image a sample and spectral information is used to
differentiate microplastic particles from nonplastic particles. In addition to
the material id, one often wants to measure the shape and size of the particles.
In a brute force technique, one could first identify every spectrum in the map,
then use thresholding and image analysis to measure the particles. However, more
often than not, particles are well separated on the image surface and background
spectra is quite different from particle spectra and therefore we can use
thresholding a priori to identify and measure the particles, then pass an
exemplary spectrum for each particle to the identification routine. It is
important to note here that this is at the bleeding edge of theory and technique
so we may be updating these functions in the near future.

## Brute Force

```{r, 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")
```

## A Priori Particle Thresholding

```{r, 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)

```

# References

Almond J, Sugumaar P, Wenzel MN, Hill G, Wallis C (2020). “Determination of the
Carbonyl Index of Polyethylene and Polypropylene Using Specified Area Under Band
Methodology with ATR-FTIR Spectroscopy.” *e-Polymers*, **20**, 369–381. doi:
[10.1515/epoly-2020-0041](https://doi.org/10.1515/epoly-2020-0041).

Chabuka BK, Kalivas JH (2020). “Application of a Hybrid Fusion Classification
Process for Identification of Microplastics Based on Fourier Transform Infrared
Spectroscopy.” *Applied Spectroscopy*, **74**(9), 1167–1183. doi:
[10.1177/0003702820923993](https://doi.org/10.1177/0003702820923993).

Cowger W, Gray A, Christiansen SH, De Frond H, Deshpande AD, Hemabessiere L, Lee
E, Mill L, et al. (2020). “Critical Review of Processing and Classification
Techniques for Images and Spectra in Microplastic Research.” *Applied
Spectroscopy*, **74**(9), 989–1010. doi:
[10.1177/0003702820929064](https://doi.org/10.1177/0003702820929064).

Cowger W, Steinmetz Z, Gray A, Munno K, Lynch J, Hapich H, Primpke S,
De Frond H, Rochman C, Herodotou O (2021). “Microplastic Spectral Classification
Needs an Open Source Community: Open Specy to the Rescue!”
*Analytical Chemistry*, **93**(21), 7543–7548. doi:
[10.1021/acs.analchem.1c00123](https://doi.org/10.1021/acs.analchem.1c00123).

Liland KH (2015). “4S Peak Filling—Baseline Estimation by Iterative Mean
Suppression.” *MethodsX*, **2**, 135–140. doi:
[10.1016/j.mex.2015.02.009](https://doi.org/10.1016/j.mex.2015.02.009).

Primpke S, Wirth M, Lorenz C, Gerdts G (2018). “Reference Database Design for
the Automated Analysis of Microplastic Samples Based on Fourier Transform
Infrared (FTIR) Spectroscopy.” *Analytical and Bioanalytical Chemistry*,
**410**(21), 5131–5141. doi: 
[10.1007/s00216-018-1156-x](https://doi.org/10.1007/s00216-018-1156-x).

Renner G, Schmidt TC, Schram J (2017). “A New Chemometric Approach for Automatic
Identification of Microplastics from Environmental Compartments Based on FT-IR
Spectroscopy.” *Analytical Chemistry*, **89**(22), 12045–12053. doi: 
[10.1021/acs.analchem.7b02472](https://doi.org/10.1021/acs.analchem.7b02472).

Savitzky A, Golay MJ (1964). “Smoothing and Differentiation of Data by
Simplified Least Squares Procedures.” *Analytical Chemistry*, **36**(8),
1627–1639.

Zhao J, Lui H, McLean DI, Zeng H (2007). “Automated Autofluorescence Background
Subtraction Algorithm for Biomedical Raman Spectroscopy.”
*Applied Spectroscopy*, **61**(11), 1225–1232. doi:
[10.1366/000370207782597003](https://doi.org/10.1366/000370207782597003).