Scatter Correction: MSC and EMSC¶

This notebook demonstrates multiplicative scatter correction (MSC) and extended multiplicative scatter correction (EMSC) for NIR-type spectra.

In [ ]:

import matplotlib.pyplot as plt
import numpy as np

from spectrakit import normalize_snv, scatter_emsc, scatter_msc
from spectrakit.plot import plot_spectrum

import matplotlib.pyplot as plt import numpy as np from spectrakit import normalize_snv, scatter_emsc, scatter_msc from spectrakit.plot import plot_spectrum

Generate Synthetic NIR-like Spectra¶

Simulate spectra with multiplicative scatter effects (common in diffuse reflectance).

In [ ]:

rng = np.random.default_rng(42)
wavenumbers = np.linspace(4000, 10000, 500)


# True underlying spectrum
def gaussian(x, c, a, s):
    return a * np.exp(-((x - c) ** 2) / (2 * s**2))


true_spectrum = (
    gaussian(wavenumbers, 5200, 0.8, 200)
    + gaussian(wavenumbers, 6900, 0.6, 300)
    + gaussian(wavenumbers, 8400, 0.4, 250)
    + 0.1
)

# Simulate 20 spectra with multiplicative scatter
n_spectra = 20
spectra = np.zeros((n_spectra, len(wavenumbers)))
for i in range(n_spectra):
    scale = rng.uniform(0.7, 1.3)
    offset = rng.uniform(-0.2, 0.2)
    noise = rng.normal(0, 0.01, len(wavenumbers))
    spectra[i] = scale * true_spectrum + offset + noise

plt.figure(figsize=(10, 4))
plot_spectrum(spectra, wavenumbers, title="Raw Spectra with Scatter Effects", invert_x=False)
plt.show()

rng = np.random.default_rng(42) wavenumbers = np.linspace(4000, 10000, 500) # True underlying spectrum def gaussian(x, c, a, s): return a * np.exp(-((x - c) ** 2) / (2 * s**2)) true_spectrum = ( gaussian(wavenumbers, 5200, 0.8, 200) + gaussian(wavenumbers, 6900, 0.6, 300) + gaussian(wavenumbers, 8400, 0.4, 250) + 0.1 ) # Simulate 20 spectra with multiplicative scatter n_spectra = 20 spectra = np.zeros((n_spectra, len(wavenumbers))) for i in range(n_spectra): scale = rng.uniform(0.7, 1.3) offset = rng.uniform(-0.2, 0.2) noise = rng.normal(0, 0.01, len(wavenumbers)) spectra[i] = scale * true_spectrum + offset + noise plt.figure(figsize=(10, 4)) plot_spectrum(spectra, wavenumbers, title="Raw Spectra with Scatter Effects", invert_x=False) plt.show()

MSC Correction¶

In [ ]:

msc_corrected = scatter_msc(spectra)

fig, axes = plt.subplots(1, 2, figsize=(14, 4))
plot_spectrum(spectra, wavenumbers, ax=axes[0], title="Before MSC", invert_x=False)
plot_spectrum(msc_corrected, wavenumbers, ax=axes[1], title="After MSC", invert_x=False)
plt.tight_layout()
plt.show()

msc_corrected = scatter_msc(spectra) fig, axes = plt.subplots(1, 2, figsize=(14, 4)) plot_spectrum(spectra, wavenumbers, ax=axes[0], title="Before MSC", invert_x=False) plot_spectrum(msc_corrected, wavenumbers, ax=axes[1], title="After MSC", invert_x=False) plt.tight_layout() plt.show()

EMSC Correction¶

EMSC adds polynomial baseline terms for wavelength-dependent scatter.

In [ ]:

emsc_corrected = scatter_emsc(spectra, poly_order=2)

fig, axes = plt.subplots(1, 2, figsize=(14, 4))
plot_spectrum(spectra, wavenumbers, ax=axes[0], title="Before EMSC", invert_x=False)
plot_spectrum(
    emsc_corrected, wavenumbers, ax=axes[1], title="After EMSC (poly_order=2)", invert_x=False
)
plt.tight_layout()
plt.show()

emsc_corrected = scatter_emsc(spectra, poly_order=2) fig, axes = plt.subplots(1, 2, figsize=(14, 4)) plot_spectrum(spectra, wavenumbers, ax=axes[0], title="Before EMSC", invert_x=False) plot_spectrum( emsc_corrected, wavenumbers, ax=axes[1], title="After EMSC (poly_order=2)", invert_x=False ) plt.tight_layout() plt.show()

MSC vs EMSC vs SNV¶

Compare scatter correction approaches.

In [ ]:

snv_corrected = normalize_snv(spectra)

fig, axes = plt.subplots(1, 3, figsize=(16, 4))

plot_spectrum(msc_corrected, wavenumbers, ax=axes[0], title="MSC", invert_x=False)
plot_spectrum(emsc_corrected, wavenumbers, ax=axes[1], title="EMSC", invert_x=False)
plot_spectrum(snv_corrected, wavenumbers, ax=axes[2], title="SNV", invert_x=False)

plt.suptitle("Scatter Correction Comparison", fontsize=14)
plt.tight_layout()
plt.show()

# Variance across spectra (lower = better correction)
print("Variance across spectra (mean across wavelengths):")
print(f"  Raw:  {spectra.var(axis=0).mean():.6f}")
print(f"  MSC:  {msc_corrected.var(axis=0).mean():.6f}")
print(f"  EMSC: {emsc_corrected.var(axis=0).mean():.6f}")
print(f"  SNV:  {snv_corrected.var(axis=0).mean():.6f}")

snv_corrected = normalize_snv(spectra) fig, axes = plt.subplots(1, 3, figsize=(16, 4)) plot_spectrum(msc_corrected, wavenumbers, ax=axes[0], title="MSC", invert_x=False) plot_spectrum(emsc_corrected, wavenumbers, ax=axes[1], title="EMSC", invert_x=False) plot_spectrum(snv_corrected, wavenumbers, ax=axes[2], title="SNV", invert_x=False) plt.suptitle("Scatter Correction Comparison", fontsize=14) plt.tight_layout() plt.show() # Variance across spectra (lower = better correction) print("Variance across spectra (mean across wavelengths):") print(f" Raw: {spectra.var(axis=0).mean():.6f}") print(f" MSC: {msc_corrected.var(axis=0).mean():.6f}") print(f" EMSC: {emsc_corrected.var(axis=0).mean():.6f}") print(f" SNV: {snv_corrected.var(axis=0).mean():.6f}")