Merge pull request 'dev/Memory_pooling_adjustments' (#1) from dev/Memory_pooling_adjustments into main

Reviewed-on: #1

.github/workflows/ci.yml

@@ -0,0 +1,32 @@
+name: CI
+
+on:
+  push:
+    branches:
+      - main
+  pull_request:
+
+jobs:
+  test:
+    name: Julia ${{ matrix.version }} - ${{ matrix.os }} - ${{ matrix.arch }} - ${{ github.event_name }}
+    runs-on: ${{ matrix.os }}
+    strategy:
+      fail-fast: false
+      matrix:
+        version:
+          - '1.10'
+        os:
+          - ubuntu-latest
+          - windows-latest
+          - macOS-latest
+        arch:
+          - x64
+    steps:
+      - uses: actions/checkout@v4
+      - uses: julia-actions/setup-julia@v2
+        with:
+          version: ${{ matrix.version }}
+          arch: ${{ matrix.arch }}
+      - uses: julia-actions/cache@v2
+      - uses: julia-actions/julia-buildpkg@v1
+      - uses: julia-actions/julia-runtest@v1

.gitignore

@@ -7,4 +7,6 @@ log/*
 R_original_scripts/
 test/results/
 test/binarization_results/
-Manifest.toml
+Manifest.toml
+MSI_sysimage.dll
+MSI_sysimage.so

Dockerfile.headless

@@ -0,0 +1,31 @@
+# Dockerfile.headless
+# This Dockerfile creates a minimal, high-speed container solely intended
+# for running JuliaMSI scripts via the pre-compiled headless sysimage.
+# It explicitly excludes the Genie UI and graphical routing to optimize size and start speed.
+
+FROM julia:1.11
+
+# System dependencies for core MSI math libraries 
+RUN apt-get update && apt-get install -y \
+    build-essential \
+    zlib1g-dev \
+    && rm -rf /var/lib/apt/lists/*
+
+WORKDIR /app
+
+# Only copy essential library files (no UI or app.jl required)
+COPY Project.toml Manifest.toml ./
+COPY src/ ./src/
+COPY build_sysimage.jl precompile_script.jl ./
+COPY test/ ./test/
+
+# Instantiate the environment
+RUN julia --project=. -e 'import Pkg; Pkg.instantiate(); Pkg.precompile()'
+
+# Build the headless sysimage
+# This command strips out Genie and plot libraries, creating a pure data-engine .so
+RUN julia --project=. build_sysimage.jl --headless
+
+# The resulting image is heavily optimized for zero-JIT execution.
+# Default command simply enters a fast REPL.
+CMD ["julia", "--project=.", "--threads", "auto", "-J", "sys_msi_headless.so"]

JuliaMSI_MainGuide.md

@@ -0,0 +1,217 @@
+# JuliaMSI Platform User Guide
+## Core Workflow and Interactive Analysis
+
+---
+
+### Table of Contents
+1. [Introduction](#introduction)
+2. [Loading a Dataset](#loading-a-dataset)
+3. [Converting mzML to imzML](#converting-mzml-to-imzml)
+4. [Batch Processing](#batch-processing)
+5. [Spectral Visualization](#spectral-visualization)
+6. [Imaging Slice Generation](#imaging-slice-generation)
+7. [Topography and Surface Plots](#topography-and-surface-plots)
+8. [Small Features and Tools](#small-features-and-tools)
+    * [Metadata Viewing](#metadata-viewing)
+    * [Optical Image Overlay](#optical-image-overlay)
+    * [Comparison of Datasets](#comparison-of-datasets)
+    * [Registry and Safe Processing](#registry-and-safe-processing)
+9. [Mask Generation and Usage (mask.jl)](#mask-generation-and-usage-maskjl)
+10. [Preprocessing](#preprocessing)
+11. [Practical Tips for Analysts](#practical-tips-for-analysts)
+    * [Julia Version Compatibility](#julia-version-compatibility)
+    * [First-Time Startup](#first-time-startup)
+    * [Accessing the Web Interface](#accessing-the-web-interface)
+    * [Memory Management](#memory-management)
+    * [Registry and File Locations](#registry-and-file-locations)
+    * [Mask Handling](#mask-handling)
+12. [Troubleshooting Common Issues](#troubleshooting-common-issues)
+13. [Citation and Contact](#citation-and-contact)
+
+---
+
+## Introduction
+This document serves as an entry point for analysts using the **JuliaMSI** platform. It focuses on the user-facing modules: `app.jl` (the main application entry), `julia_imzml_visual.jl` (visualization and heavy procedures), and `mask.jl` (mask generation and filtering). The goal is to explain why each step requires specific inputs and what those inputs mean scientifically.
+
+**JuliaMSI** is a user-friendly Graphical User Interface (GUI) developed in Julia for MSI data analyses. Our framework allows for rapid **imzML** MSI data analysis and **mzML** spectra plot loading, significantly improving analytical workflows. High-resolution instruments produce massive datasets, requiring efficient memory management and parallel processing—reasons why JuliaMSI is written in Julia, which offers high performance.
+
+### Core Objectives:
+*   **Simplification**: Reducing command-line dependency through intuitive visualizations and easy parameter input.
+*   **Accessibility**: Breaking the proprietary software barrier across fields like plant sciences, biomedical research, and environmental studies.
+*   **Interactive Exploration**: Using the `Plotly.jl` package for dynamic visualizations. Users can explore multiple interpretations of images, select specific coordinates for spectra analysis, or pick an *m/z* value directly from the displayed spectrum.
+
+### Key Capabilities:
+*   **Smart data handling**: Processes large files in manageable sections without overwhelming memory.
+*   **Batch processing**: Analyzes multiple samples simultaneously for high-throughput studies.
+*   **Interactive exploration**: Web interface enables real-time data visualization and analysis.
+
+---
+
+## Loading a Dataset
+**Action**: Click the **Search** button (`@onbutton btnSearch`)
+
+Load an **imzML** or **mzML** file into the workspace, extract metadata, and prepare the dataset for visualization and processing.
+
+The platform checks whether a pre-processed version of the same file exists in the registry (**fast load**). If not, it performs a **full load** using `OpenMSIData`. The GUI handles multiple file uploads for MSI slice generation; however, for spectra and preprocessing, it will only use the last selected file.
+
+---
+
+## Converting mzML to imzML
+**Action**: Click **Convert/Process** (`@onbutton convert_process`)
+
+This tool converts an **mzML** file (from a non-imaging mass spectrometry run) into an **imzML** file suitable for imaging analysis. This is necessary when spectra were acquired in a raster pattern but saved as a continuous **mzML** file without spatial coordinates.
+
+### Inputs Required:
+*   **mzML file**: The raw mass spectrometry data file.
+*   **Sync file**: A text file that maps each spectrum in the **mzML** to a specific pixel coordinate (*x, y*). The sync file typically contains one line per spectrum with the format `x y` or similar.
+
+> **Why these inputs?**: **mzML** files lack spatial metadata; they contain a list of spectra in acquisition order. Without coordinate information, the data cannot be reconstructed as an image. **Sync files** provide the missing spatial mapping, often generated by instrument software from a motor stage log. The conversion process embeds these coordinates into the **imzML** structure.
+
+---
+
+## Batch Processing
+**Action**: Click **Main Process** (`@onbutton mainProcess`)
+
+This is the core image generation routine. It takes a list of comma-separated *m/z* values and produces ion images (and optionally **TrIQ Quantized** images) for each file loaded.
+
+### Required Inputs:
+*   **m/z values**: A comma-separated list of *m/z* values to visualize spatial distribution. These are the core of MSI: each image maps the abundance of a particular ion across the sample. Without a target, there is no image.
+*   **Tolerance**: The mass window around each target ion. Necessary because no mass spectrometer measures masses with infinite precision. The signal from a given ion is spread over a small *m/z* range; integrating over that range (by summing intensities of all data points within ±tolerance) captures the true signal.
+*   **Color level**: Number of discrete color levels in the output BMP image. A low number (e.g., 16) gives a coarse view; a high number (256) shows subtle intensity variations but can be noisy.
+*   **TrIQ**: A toggle to enable **Threshold Intensity Quantization**. This produces images with enhanced contrast and reduced influence of extreme outliers, making spatial patterns more visible. It is a contrast normalization and quantization technique.
+*   **TrIQ Threshold**: Probability threshold for outlier removal in TrIQ.
+*   **Median filter**: Replaces each pixel’s intensity with the median of its neighbors, effectively removing isolated "hot pixels" while preserving sharp boundaries.
+*   **Mask**: Toggle to focus on biologically relevant regions (e.g., tissue vs. background), avoiding wasted computational resources on irrelevant pixels.
+
+> **Note**: The mask path is determined inside the JSON file, which can be manually altered, or using the specific section in `mask.jl`.
+
+---
+
+## Imaging Slice Generation (Part of Main Process)
+In MSI, each pixel contains a full mass spectrum. To create an ion image (or "slice"), we must decide which data points in that spectrum belong to the ion of interest.
+
+*   **Tolerance**: Defines the *m/z* window. Too narrow may miss part of the peak; too wide includes noise or interfering peaks. Optimal tolerance depends on instrument resolving power and calibration stability.
+*   **Colorbar Generation**: Each image is accompanied by a colorbar PNG showing the mapping between intensity and colors. The bounds are determined during quantization (TrIQ or min-max scaling) and passed to `generate_colorbar_image`. This ensures consistent color interpretation with labels indicating the intensity range and scaling factors (e.g., "Intensity ×10³").
+
+---
+
+## Spectral Visualization
+**Functions**: `createMeanPlot`, `createSumPlot`, `createXYPlot`, `createNSpectrumPlot`
+
+Displays spectra from the loaded dataset, optionally filtered by a mask. This generates a Plotly layout plot for quality control, identification, and *m/z* selection.
+
+### Inputs & Features:
+*   **Dataset**: The last **imzML** / **mzML** file loaded in the GUI.
+*   **Mask**: Toggle for mask filtering. The spectrum is computed only from pixels inside the mask.
+*   **User click**: If a user clicks on an already generated **slice** (2D image), two vertical lines appear and values are filled (one for *x*-axis, one for *y*-axis) to generate spectrum plot indicating the position of the click in the image. This allows the creation of a spectrum plot correlating to specific spatial selection.
+*   **NSpectrumPlot**: Allows the user to select specific spectra based on the **ID** (the index of the spectra in the dataset).
+
+### Spectral Types:
+*   **Mean spectrum**: The average of all pixel spectra. Gives an overview of ions present in the entire sample and their relative abundances. Used to identify peaks for further imaging.
+*   **Sum spectrum**: Total ion current summed across all pixels. Emphasizes ions that are abundant across many pixels.
+*   **Per-pixel spectra (XY and N)**: Reveal the distribution of ions in a specific location, useful for comparing tissue regions (e.g., tumor vs. normal).
+
+> **Why these inputs?**: Spectra are the raw data of MSI. Visualizing them is the first step: checking for unexpected peaks, assessing signal-to-noise ratio, and verifying mass calibration. Applying a mask (e.g., only tumor region) allows researchers to extract region-specific molecular signatures.
+
+---
+
+## Topography and Surface Plots
+**Functions**: `image3dPlot`, `imageCPlot`, `triq3dPlot`, `triqCPlot`
+
+Convert a 2D ion image into a 3D surface plot where the height (*z*-axis) represents ion intensity.
+*   **Inputs**: Current ion image and optional Mask.
+*   **Why these inputs?**: 3D surfaces and 2D Contour plots make intensity variations more tangible, especially for features with large dynamic range. They are useful for understanding the topography of ion distributions (e.g., a lipid enriched in a specific tissue layer). The 3D plot is interactive, allowing for rotation and zoom.
+
+---
+
+## Small Features and Tools
+
+### Metadata Viewing (`showMetadataBtn`)
+Display a summary of the dataset’s metadata: instrument settings, acquisition parameters, image dimensions, number of pixels, *m/z* range, etc.
+
+### Optical Image Overlay (`btnOptical`)
+Overlay a histological or photographic image onto the MSI image to correlate molecular distributions with tissue anatomy.
+*   **Inputs**: Current ion image, Transparency slider (`imgTrans`), and Optical image file (PNG, BMP, JPG, JPEG).
+*   **Output**: The main image tab updates to show the MSI image (or TrIQ image) blended with the optical background.
+
+### Comparison of Datasets
+Functionality to compare two datasets (e.g., control vs. treated). This involves side-by-side visualization or ratio image generation to highlight differences in molecular abundance. It allows visualization of slices with or without TrIQ filtering.
+
+### Registry and Safe Processing
+Batch processing is orchestrated by `process_file_safely`. The registry (`registry.json`) stores critical info: source path, image dimensions, summary statistics, and mask path. This file can be edited manually following its generated structure.
+
+---
+
+## Mask Generation and Usage (`mask.jl`)
+Create binary masks that define **Regions of Interest (ROI)** in the sample. Masks filter spectra, preprocessing, and images to include only pixels inside the ROI.
+
+### User-Defined Inputs:
+*   **Drawing tools**: Draw freehand shapes on the ion image or optical overlay.
+*   **Otsu Scale**: Automatically generates a mask based on the Otsu thresholding algorithm with the slider percentage.
+*   **Noise Size**: Removes noise from the mask based on the slider percentage.
+*   **Hole Size**: Removes holes in the mask based on the slider percentage.
+*   **Smoothing**: Smooths the mask based on the slider percentage.
+*   **Thresholding**: Set intensity thresholds to automatically segment regions.
+*   **Import**: Load a pre-existing mask image (e.g., from a segmentation tool).
+
+### Outputs:
+*   A binary mask matrix (same dimensions as the MSI image) saved in the public/css/masks folder as a `.png` file with the name of the dataset.
+*   The mask is registered in the dataset’s metadata entry for automatic application in subsequent analyses.
+
+> **Why these inputs?**: Often only a portion of the imaged area contains tissue; the rest could be background. Masking removes background pixels, improving statistical power and avoiding artifacts. Manual delineation based on optical images is standard practice in histology-directed MSI. It can also work for a selection of an area of interest (e.g. only the tumor region) in samples with multiple areas of interest.
+
+---
+
+## Preprocessing
+To understand the deep structure of the preprocessing pipeline, please refer to the dedicated guide: **[JuliaMSI Preprocessing Guide](JuliaMSI_Preprocessing.md)**
+
+This document describes the rationale, methods (TIC/PQN, SNIP, Wavelets), and implementation choices in detail.
+
+---
+
+## Practical Tips for Analysts
+
+### Julia Version Compatibility
+*   The platform is tested and compatible with **Julia 1.11**.
+*   **Recommendation**: Use Julia 1.11 to ensure all dependencies work correctly. Use **Juliaup** to manage multiple versions easily. Downgrading your system Julia is unnecessary; you can keep multiple versions and launch the platform with the specific executable path.
+
+### First-Time Startup
+*   The script will detect a missing `Manifest.toml` and automatically instantiate packages. This may take a few minutes. Subsequent startups will be much faster.
+*   **Caution**: Do not close the terminal/console while the server is running; it hosts the web interface.
+
+### Accessing the Web Interface
+*   The platform opens your browser to `http://127.0.0.1:1481`.
+*   If port 1481 is in use, change this segment in the startup script: `up(host="127.0.0.1", port=1481)`.
+
+### Memory Management
+*   MSI datasets can be large (gigabytes). JuliaMSI includes automatic **garbage collection** and, on Linux, calls `malloc_trim` to return freed memory to the OS.
+*   **OOM Troubleshooting**: Process fewer *m/z* values at once, use a mask to exclude background, or increase system swap space.
+
+### Registry and File Locations
+*   `registry.json` stores paths to your datasets. If you move original **imzML** files, update the paths in the JSON or delete the registry file to let the platform rebuild it.
+*   Processed images and colorbars are saved in the `public/` folder, organized by dataset name.
+
+### Mask Handling
+*   Ensure manually created masks have the **exact same dimensions** (width × height) as the MSI image.
+*   Masks are automatically applied during batch processing if the "Mask" toggle is enabled. If you move the mask file, you must update its path in the registry.
+
+---
+
+## Troubleshooting Common Issues
+
+*   **"Error loading active file"**: Check that the file is a valid **imzML** or **mzML** and that the corresponding `.ibd` file is in the same directory.
+*   **Images appear all black**: This can happen if the tolerance is too narrow (no signal captured) or if the mask is incorrectly applied. Verify *m/z* and tolerance.
+*   **Slow performance**: Reduce the number of *m/z* values or disable the median filter. Ensure you have enough RAM.
+*   **Browser interface not responding**: Refresh the page. Check the console for errors; you may need to restart the Julia process.
+
+---
+
+## Citation and Contact
+Any revision to the code, bug, or error should be reported to: **julian.sierrag@icloud.com**
+
+*Note: This guide is based on collective knowledge of these algorithms; please perform your own research and cite properly. To cite the JuliaMSI platform, please use:*
+
+> José Julián Sierra-Álvarez, Martín Orlando Camargo-Escalante, Carlos Daniel Sierra-Álvarez, Carmelo Hernández-Caricio, Juan Francisco Moreno-Luna, Isabel Buendía-Corona, Robert Winkler,
+> **JuliaMSI: A high-performance graphical platform for mass spectrometry imaging data analysis**,
+> *Analytica Chimica Acta*, Volume 1377, 2025, 344613, ISSN 0003-2670,
+> [https://doi.org/10.1016/j.aca.2025.344613](https://doi.org/10.1016/j.aca.2025.344613)

JuliaMSI_Preprocessing.md

@@ -0,0 +1,354 @@
+# A Comprehensive Guide to Mass Spectrometry Imaging Preprocessing
+## Rationale, Methods, and Implementation
+
+---
+
+### Table of Contents
+1. [Introduction](#introduction)
+2. [Profile vs. Centroid Mode Data](#profile-vs-centroid-mode-data)
+3. [General Principles of MSI Preprocessing](#general-principles-of-msi-preprocessing)
+4. [Step-by-Step Breakdown of Preprocessing Steps](#step-by-step-breakdown-of-preprocessing-steps)
+    * [Validation and Quality Control](#validation-and-quality-control)
+    * [Intensity Transformation (Stabilization)](#intensity-transformation-stabilization)
+    * [Smoothing](#smoothing)
+    * [Baseline Correction](#baseline-correction)
+    * [Normalization](#normalization)
+    * [Peak Detection](#peak-detection)
+    * [Peak Selection (Filtering)](#peak-selection-filtering)
+    * [Calibration](#calibration)
+    * [Peak Alignment](#peak-alignment)
+    * [Peak Binning and Feature Matrix Generation](#peak-binning-and-feature-matrix-generation)
+5. [Example Walkthrough: From Raw Spectra to Feature Matrix](#example-walkthrough-from-raw-spectra-to-feature-matrix)
+    * [Phase 1 – Instrument & Data Characteristics](#phase-1--instrument--data-characteristics)
+    * [Phase 2 – Noise & Signal Quality](#phase-2--noise--signal-quality)
+    * [Phase 3 – Mass Accuracy Analysis](#phase-3--mass-accuracy-analysis)
+    * [Phase 4 – Spatial Region Analysis](#phase-4--spatial-region-analysis)
+    * [Phase 5 – Peak Characteristics](#phase-5--peak-characteristics)
+    * [Phase 6 – Preprocessing Recommendations](#phase-6--preprocessing-recommendations)
+6. [Extra Tips for the Analyst](#extra-tips-for-the-analyst)
+    * [Isotopic and Charge State Matters](#isotopic-and-charge-state-matters)
+    * [Always Visualize the Average Spectrum First](#always-visualize-the-average-spectrum-first)
+    * [Preprocessing Is Not a One-Shot Process](#preprocessing-is-not-a-one-shot-process)
+
+---
+
+## Introduction
+Mass Spectrometry Imaging (MSI) generates thousands of mass spectra, each acquired from a distinct pixel location. In the literature, before any interpretation can be drawn from these data, a series of preprocessing steps must be applied. These steps correct for instrumental artifacts, reduce noise, and align data across pixels to produce a clean, comparable feature matrix suitable for statistical analysis and image reconstruction.
+
+Inside the **JuliaMSI** code (`Preprocessing.jl` and `PreprocessingPipeline.jl`), we have a preprocessing workflow inspired by established R packages like **MALDIquant** and **Cardinal**.
+
+This document explains the rationale behind each step, the required inputs, the scientific meaning of parameters, and the context behind key algorithms such as **TIC normalization**, **PQN**, **wavelet-based peak detection**, and **LOWESS alignment**.
+
+---
+
+## Profile vs. Centroid Mode Data
+
+Mass spectrometry data can be acquired in two primary modes: **profile** and **centroid**. Understanding the difference is crucial because preprocessing steps must be adapted accordingly.
+
+*   **Profile mode**:
+    The instrument records intensity at every digitized point along the *m/z* axis, producing a continuous waveform. Peaks appear as Gaussian-like shapes. This mode retains full peak shape information, which is valuable for accurate mass measurement, resolution assessment, and advanced peak detection. However, file sizes are large and noise is present between peaks.
+    
+*   **Centroid mode**:
+    The instrument performs real-time peak detection and reduces each peak to a single centroid (*m/z*, intensity) pair, discarding the underlying profile. This dramatically reduces file size but loses peak shape details. Centroid data are essentially a list of detected peaks with no information about peak width or baseline.
+
+In **profile mode**, algorithms must operate on the full *m/z* and intensity vectors, performing tasks like smoothing, baseline correction, and peak picking on continuous data. In **centroid mode**, these steps are often unnecessary or simplified: smoothing and baseline correction are not applicable because the data are already peak lists; peak picking is replaced by filtering based on intensity thresholds. 
+
+> [!IMPORTANT]
+> The JuliaMSI code internally distinguishes this difference in some steps by acquiring the mode from the metadata on the files, but it's up to the analyst to decide skipping smoothing/baseline steps for centroid data, for example.
+
+---
+
+## General Principles of MSI Preprocessing
+
+All preprocessing steps share a common goal: **to end up with an aligned feature matrix where rows correspond to pixels (spectra) and columns correspond to *m/z* bins (features).**
+
+The steps are applied sequentially, often in a specific order to avoid compounding errors. For example, **baseline correction should precede normalization** because baseline offsets can distort total ion current calculations. Similarly, **peak alignment should follow peak detection** because alignment requires a set of common peaks across spectra.
+
+---
+
+## Step-by-Step Breakdown of Preprocessing Steps
+
+### Validation and Quality Control
+*(This step doesn't require analyst input)*
+*   **Functions in code**: `validate_spectrum`, `qc_is_empty`, `qc_is_regular`
+*   **Purpose**: Ensure that each input spectrum is well-formed and contains usable data before any processing.
+*   **Inputs required**:
+    *   `mz::AbstractVector{<:Real}`: The *m/z* values.
+    *   `intensity::AbstractVector{<:Real}`: The corresponding intensities.
+
+> **Why these inputs?**
+> Every subsequent algorithm depends on the integrity of these vectors. Checking for non-emptiness, equal length, finite values, non-negativity, and monotonic *m/z* order prevents cryptic errors later.
+
+---
+
+### Intensity Transformation (Stabilization)
+*   **Function**: `transform_intensity_core`
+*   **Purpose**: Apply a variance-stabilizing transformation to make the data more homoscedastic (constant variance across intensity range).
+*   **Inputs required**:
+    *   `intensity::AbstractVector{<:Real}`
+    *   `method::Symbol` (e.g., `:sqrt`, `:log1p`)
+
+> **Why these inputs?**
+> The transformation operates on the intensity vector only; *m/z* values remain unchanged. The choice of method depends on the data characteristics (e.g., presence of zeros).
+>
+> Stabilization uses functions like **square root**, **log**, **log2**, **log10** and **log1p** helps the analyst stabilize variance for count data, or handle quadratical variance. The analyst must know which one of these methods to use according to the data they are visualizing (e.g., by the mean spectrum plot, or a plot from a specific coordinate).
+>
+> In MSI in general, **square root stabilization** is common because ion counting statistics approximate Poisson.
+>
+> **Example**: A spectrum with intense peaks and low baseline noise: after square root, the dynamic range is compressed, making low-intensity features more comparable to high-intensity ones.
+
+---
+
+### Smoothing
+*   **Function**: `smooth_spectrum_core`, `moving_average_smooth`
+*   **Purpose**: Reduce high-frequency noise while preserving peak shapes.
+*   **Inputs required**:
+    *   `y::AbstractVector{<:Real}`: Intensity vector.
+    *   `method::Symbol` (`:savitzky_golay`, `:moving_average`)
+    *   `window::Int`: Size of the smoothing window.
+    *   `order::Int`: Polynomial order for Savitzky–Golay.
+
+> **Why these inputs?**
+> Smoothing operates on the intensity vector only. The window size must be an odd integer for symmetry; the polynomial order must be less than the window. These constraints ensure the filter behaves properly (e.g., Savitzky–Golay requires the window bigger than the order parameter).
+>
+> **Savitzky–Golay** fits a low-degree polynomial to the local window, preserving moments of the peak (height, width) better than a simple moving average, which can distort peak shapes.
+>
+> **Example**: A noisy peak with 5% random noise: after Savitzky–Golay, the peak becomes smoother while its centroid and height remain accurate.
+
+---
+
+### Baseline Correction
+*   **Function**: `apply_baseline_correction_core`, `_snip_baseline_impl`, `convex_hull_baseline`, `median_baseline`
+*   **Purpose**: Estimate and subtract the background signal (baseline) from each spectrum.
+*   **Inputs required**:
+    *   `y::AbstractVector{<:Real}`: Intensity vector.
+    *   `method::Symbol` (`:snip`, `:convex_hull`, `:median`)
+    *   `iterations::Int` (for SNIP)
+    *   `window::Int` (for median)
+
+> **Why these inputs?**
+> Baseline correction requires only the intensity vector. The method-specific parameters control the aggressiveness and shape of the baseline estimate.
+>
+> *   **SNIP (Sensitive Nonlinear Iterative Peak clipping)**: Iteratively replaces each point with the minimum of itself and the average of its neighbors. After many iterations, only the baseline remains. SNIP is the default in **MALDIquant** due to its effectiveness.
+> *   **Convex hull**: Constructs the upper convex hull of the spectrum (or lower hull for baseline) and interpolates.
+> *   **Moving median**: Estimates baseline as the local median. Simple but can be biased by dense peaks.
+>
+> **Example**: In a spectrum with a rising baseline after SNIP correction, the baseline is removed, revealing peaks sitting on a flat background.
+
+---
+
+### Normalization
+*   **Function**: `apply_normalization_core`, `tic_normalize`, `pqn_normalize`, `median_normalize`, `rms_normalize`
+*   **Purpose**: Correct for differences in total ion abundance between spectra, making intensities comparable across pixels.
+*   **Inputs required**:
+    *   `y::AbstractVector{<:Real}` (per-spectrum) or a matrix `M` (for PQN).
+    *   `method::Symbol`
+
+> **Why these inputs?**
+> Normalization operates on intensities. **PQN** requires a matrix of spectra (rows = *m/z* bins, columns = spectra) because it computes a reference spectrum across all samples.
+>
+> *   **TIC (Total Ion Current)**: Divides each intensity by the sum of all intensities in the spectrum. A simple algorithm but sensitive to a few very intense peaks (e.g., matrix peaks).
+> *   **Median normalization**: Divides by the median intensity. Robust to outliers but can be unstable if many intensities are zero.
+> *   **RMS normalization**: Divides by root-mean-square, preserving the shape of the spectrum while scaling to unit energy.
+> *   **PQN (Probabilistic Quotient Normalization)**: Calculates a reference spectrum and scales each spectrum by the median of the quotients (intensity/reference). This method is preferred when some peaks vary drastically because it downweights their influence.
+>
+> **Example**: Two spectra from different tissue regions: one has a huge lipid peak, the other has lower overall signal. TIC would suppress the lipid peak in the first spectrum, possibly losing biological information. PQN would scale both spectra more equitably.
+
+---
+
+### Peak Detection
+*   **Functions**: `detect_peaks_profile_core`, `detect_peaks_wavelet_core`, `detect_peaks_centroid_core`
+*   **Purpose**: Identify peaks (signals of interest) in each spectrum.
+*   **Inputs required**:
+    *   *m/z*, intensity vectors.
+    *   Method-specific parameters: `snr_threshold`, `half_window`, `min_peak_prominence`, `merge_peaks_tolerance`, etc.
+
+> **Why these inputs?**
+> Peak detection uses both *m/z* and intensity to locate local maxima and evaluate peak quality. Parameters control sensitivity and filtering. A peak is a local maximum that rises significantly above the local noise. Detection algorithms vary by data mode:
+>
+> *   **Profile mode (local maxima)**: Scans for points that are highest in a sliding window (`half_window`). Then filters by **SNR** (signal-to-noise ratio) and prominence (height relative to surrounding minima). Close peaks are merged if within `merge_peaks_tolerance`.
+> *   **Wavelet-based**: Uses **Continuous Wavelet Transform (CWT)** to identify peaks at multiple scales, making it robust to varying peak widths. Local maxima in the CWT magnitude indicate peak positions. Wavelet is often used for complex spectra with overlapping peaks. CWT-based peak detection is implemented in **MALDIquant** as the default method for profile data.
+> *   **Centroid mode**: Simply filters the existing peak list by SNR, as peaks are already identified.
+>
+> **Example**: A spectrum with a small peak on the shoulder of a large peak: wavelet detection may resolve both, while simple local maxima might miss the shoulder.
+
+---
+
+### Peak Selection (Filtering)
+*   **Function**: `apply_peak_selection`
+*   **Purpose**: Remove peaks that do not meet quality criteria, such as minimum SNR, acceptable width, or good Gaussian shape.
+*   **Inputs required**:
+    *   peaks list (from peak detection).
+    *   Thresholds: `min_snr`, `min_fwhm_ppm`, `max_fwhm_ppm`, `min_shape_r2`.
+
+> **Why these inputs?**
+> These thresholds reflect the analyst's confidence in peak quality. They are typically derived from instrument specifications or empirical data.
+>
+> Not all local maxima are true peaks. Noise spikes, baseline artifacts, or electronic interference can produce spurious peaks. Filtering by **SNR** removes noise. **FWHM** (full width at half maximum) in **ppm** ensures that peaks have physically reasonable widths (e.g., isotopic peaks have predictable spacing). **Shape R²** from a Gaussian fit assesses how well the peak resembles a theoretical peak; poor shape may indicate saturation or interference.
+>
+> **Example**: A peak with SNR < 3 is likely noise and is discarded. A peak with FWHM > 100 *ppm* might be background rather than a specific ion.
+
+---
+
+### Calibration
+*(This step requires that you have a list of Internal Standards inside JuliaMSI preprocessing tab filled in)*
+*   **Function**: `calibrate_spectra_core`, `find_calibration_peaks_core`
+*   **Purpose**: Correct systematic mass errors by aligning the *m/z* axis to known reference masses (internal standards).
+*   **Inputs required**:
+    *   spectra (list of `MutableSpectrum`).
+    *   `internal_standards`: dictionary of theoretical *m/z* values.
+    *   `ppm_tolerance`: matching tolerance.
+    *   `fit_order`: polynomial order for calibration curve.
+
+> **Why these inputs?**
+> Calibration requires a set of known masses that are present in the spectra. The algorithm matches detected peaks to these references, then fits a transformation (e.g., linear) to map measured to theoretical *m/z*.
+>
+> Mass spectrometers drift; this step is made to correct this drift, ensuring that the same *m/z* value corresponds to the same ion across the entire imaging dataset. **Linear or quadratic calibration** is standard. The code currently uses linear interpolation between matched peaks, which is sufficient for most MSI applications where drift is small.
+>
+> **Example**: A reference peak expected at **760.585 Da** (a lipid) is found at **760.590 Da** in a spectrum. Calibration shifts the entire axis so that the peak aligns with the theoretical value.
+
+---
+
+### Peak Alignment
+*   **Function**: `align_peaks_lowess_core`
+*   **Purpose**: Correct residual *m/z* shifts between spectra after calibration, ensuring that the same analyte peak appears at the same *m/z* across all pixels.
+*   **Inputs required**:
+    *   Reference spectrum's peak *m/z* list.
+    *   Target spectrum's peak *m/z* list.
+    *   Parameters: method (`:lowess`, `:linear`, `:ransac`), `tolerance`, `tolerance_unit`, `max_shift_ppm`.
+
+> **Why these inputs?**
+> Alignment operates on peak lists, not full spectra, to be efficient. It finds common peaks between a reference spectrum (e.g., the one with most peaks) and each target, then fits a warping function.
+>
+> Even after calibration, small nonlinear shifts may persist due to sample topography or matrix effects. Alignment brings all spectra onto a common *m/z* scale, which is essential for constructing a feature matrix.
+>
+> **LOWESS (Locally Weighted Scatterplot Smoothing)**: Fits a smooth curve through the matched peak pairs, allowing for nonlinear corrections. It is robust to outliers and does not assume a global polynomial shape. It has been used in **MALDIquant**. **RANSAC (Random Sample Consensus)** is a robust alternative. The code defaults to LOWESS because it is well-suited for smooth, continuous drift.
+>
+> **Example**: Peaks in a target spectrum are systematically shifted by **+0.01 Da** at low *m/z* and **+0.005 Da** at high *m/z*. LOWESS fits a curve that corrects each peak individually, producing aligned *m/z* values.
+
+---
+
+### Peak Binning and Feature Matrix Generation
+*   **Function**: `bin_peaks_core`
+*   **Purpose**: Group peaks from all spectra into common *m/z* bins, producing a final feature matrix.
+*   **Inputs required**:
+    *   spectra (with detected peaks).
+    *   `params::PeakBinning`: method (`:adaptive` or `:uniform`), `tolerance`, `frequency_threshold`, etc.
+
+> **Why these inputs?**
+> Binning requires the list of peaks from each spectrum. The method determines how bins are defined.
+>
+> *   **Adaptive binning**: Clusters peaks based on their *m/z* proximity (using tolerance). Bins are created where peaks are dense, and bins without enough peaks are discarded.
+> *   **Uniform binning**: Divides the total *m/z* range into a fixed number of equally spaced bins. 
+>
+> After alignment, peaks from different spectra still have slightly different *m/z* values due to residual noise. Binning assigns them to discrete bins, creating a matrix suitable for multivariate analysis. The **frequency threshold** removes bins that appear in too few spectra (likely noise).
+>
+> **Example**: Peaks around *m/z* **760.58** from 100 spectra are grouped into one bin. The intensity for each spectrum in that bin is the maximum intensity of peaks assigned to it.
+
+---
+
+## Example Walkthrough: From Raw Spectra to Feature Matrix
+### with Intelligent Parameter Recommendations
+
+To illustrate the entire pipeline, consider a mass spectrometry imaging dataset of a mouse brain section acquired in **profile mode**. The analyst has a list of internal standards for lipids (e.g., known phospholipids at *m/z* **760.585**, **798.540**, **834.528**) to aid in mass calibration.
+
+Before running the full preprocessing, the **JuliaMSI** code executes the intelligent pre-analysis pipeline (`main_precalculation` which is a wrapper for `run_preprocessing_analysis`). This function samples a representative subset of spectra (e.g., 100 pixels) and performs a battery of diagnostic tests to suggest optimal parameters.
+
+#### The pre-analysis yields a report:
+
+*   **Phase 1 – Instrument & Data Characteristics**:
+    *   Acquisition mode: **profile** (continuous waveforms)
+    *   *m/z* axis type: **regular** (constant step of ~0.01 Da)
+    *   Dynamic range: ~4 orders of magnitude
+*   **Phase 2 – Noise & Signal Quality**:
+    *   Estimated noise level (median absolute deviation): 45 counts
+    *   Suggested SNR threshold: **3.0**
+    *   TIC coefficient of variation: **35%** → indicates moderate heterogeneity, so robust normalization (e.g., **PQN**) may be beneficial.
+*   **Phase 3 – Mass Accuracy Analysis** (using the provided lipid standards):
+    *   Mean PPM error: **12.3 ppm**
+    *   Standard deviation: **8.1 ppm**
+    *   Suggested bin tolerance (mean + 3σ): **36.6 ppm** (capped between 10 and 100 ppm).
+*   **Phase 4 – Spatial Region Analysis**:
+    *   *(if region masks were supplied) – skipped in this example.*
+*   **Phase 5 – Peak Characteristics**:
+    *   Mean FWHM: **22.5 ppm** (converted from Δm at *m/z* 500)
+    *   Mean Gaussian R²: **0.85** (peaks are well-shaped)
+    *   Mean peaks per spectrum: **187**
+*   **Phase 6 – Preprocessing Recommendations**:
+    *   The system synthesizes all information into a dictionary of suggested parameters for each step. Parameters that cannot be determined automatically (e.g., frequency threshold for binning, minimum matched peaks) are left as `nothing`; the analyst must supply them based on experimental goals.
+
+> [!NOTE]
+> If you have internal standards, remember to re-run the preprocessing recommendations after you have loaded them. This is only for precalculation and is not a required step for analysts that know about their data.
+
+#### Execution of the Pipeline:
+The analyst can now execute the pipeline, possibly overriding a few parameters (e.g., increasing SNR threshold to 4 for more stringent peak picking).
+
+1.  **Validation**: All spectra pass QC; one empty spectrum is flagged and removed.
+2.  **Intensity transformation (sqrt)**: Stabilizes variance, making low-intensity features more comparable.
+3.  **Baseline correction (SNIP, 100 iterations)**: Removes background using the recommended window of 200 points.
+4.  **Smoothing (Savitzky–Golay, window 9, order 2)**: Reduces noise while preserving peak shape.
+5.  **Peak detection (profile mode)**: Detects peaks with SNR ≥ 3, prominence ≥ 0.002, and uses the recommended half-window of 5 points. Peaks closer than 18 *ppm* are merged.
+6.  **Peak selection**: Filters peaks: keeps only those with SNR ≥ 3, FWHM between 11 and 55 *ppm*, and shape R² ≥ 0.68.
+7.  **Calibration**: Using the internal standards, each spectrum’s *m/z* axis is linearly corrected based on matched peaks within 36.6 *ppm* tolerance.
+8.  **Peak alignment (LOWESS)**: Aligns all spectra to a reference spectrum (the one with most peaks) with a span of 0.6, using the same 36.6 *ppm* tolerance. This corrects residual nonlinear drift.
+9.  **Normalization (RMS)**: Divides each spectrum by its root-mean-square intensity, chosen because TIC variation was high.
+10. **Peak binning (adaptive)**: Groups peaks within 36.6 *ppm* into bins. Bins containing fewer than 3 peaks or appearing in <20% of spectra are discarded. The final bin centers are intensity-weighted averages.
+11. **Feature matrix**: A **999 × 425 matrix** is produced, where rows are pixels and columns are *m/z* bins. This matrix is ready for PCA, segmentation, or statistical testing.
+
+#### How the Analyst Uses the Recommendations:
+*   The suggested values are data-driven, but the analyst can override them based on domain knowledge. Always verify; experimenting is encouraged as data integrity is kept.
+*   Parameters like `frequency_threshold` in peak selection and binning require the analyst’s judgment. A threshold of **0.2** means that a bin must contain a peak in at least 20% of all pixels to be retained. This reduces noise but may discard rare but real signals.
+*   This is an iterative refinement; the analyst may rerun the analysis with adjusted parameters (e.g., a lower bin tolerance).
+*   JuliaMSI aims for a tailored preprocessing workflow that respects both the instrument’s characteristics and the biological question.
+
+---
+
+## Extra Tips for the Analyst
+
+The preprocessing pipeline is powerful, but its success ultimately depends on the analyst’s informed choices. Here are practical tips to guide you through the process, especially when working with internal standards and interpreting intermediate results.
+
+> [!WARNING]
+> This is a really resource-intensive pipeline dependent on the sample, so run your precautions.
+
+### Isotopic and Charge State Matters
+When providing a dictionary of reference masses for calibration, ensure that the theoretical *m/z* values correspond to the same charge state and isotopic composition as the peaks you expect to see in your data.
+*   **Example**: If your instrument detects lipids mainly as **[M+H]⁺** ions, the reference mass must be the monoisotopic mass plus a proton. Using the neutral mass will introduce a systematic offset of **1.0078 Da**, leading to poor calibration. (Make the same adjustment for negative ion mode).
+*   If isotopic peaks are prominent, you may want to include multiple isotopes (e.g., **M**, **M+1**, **M+2**) as separate references, but be aware that their relative intensities vary. The calibration algorithm will match the closest detected peak. This is common for MSI preprocessing.
+*   The pre-analysis suggests a bin tolerance based on observed mass error. A tolerance that is too large may match the wrong peak (e.g., an isobaric interference); a tolerance that is too small may miss valid matches.
+
+### Always Visualize the Average Spectrum First
+Before running any preprocessing, generate a sum or average spectrum. This single plot tells you:
+1.  **Baseline shape**: If it's flat, rising, or undulating, it can guide the baseline correction method (**SNIP** works well for most shapes; **convex hull** is faster but may fail on curvy baselines).
+2.  **Noise level**: By zooming in.
+3.  ***m/z* Scatter**: How scattered are the spectra *m/z*.
+4.  **Peak density & Intensity range**: Influence the peak detection thresholds (SNR, prominence) and binning strategy.
+
+*   If the average spectrum shows a **rising baseline** from low to high *m/z*, SNIP with ~100 iterations will handle it well. If extremely wavy, increase iterations or try a different method.
+*   If running from the console (`run_preprocessing.jl`), the output of `run_preprocessing_analysis` provides:
+    *   **TIC coefficient of variation (CV)**: A CV > 30% suggests substantial pixel-to-pixel variation, favoring **RMS** or **PQM** over TIC.
+    *   **Mean PPM error**: If >20 *ppm*, calibration is essential. If <5 *ppm* and drift is linear, you might skip alignment.
+    *   **FWHM and Gaussian R²**: Indicate peak quality. If R² < 0.6, peaks may be poorly shaped; consider increasing the `min_peak_shape_r2` threshold.
+*   **Acquisition mode**: For **centroid data**, never apply smoothing or baseline correction because those steps assume continuous waveforms.
+    *   The pipeline will set `:method = :centroid` in peak picking and use permissive filters.
+    *   Use a low SNR threshold (e.g., 2.0).
+    *   Disable shape-based filtering (`min_peak_shape_r2 = 0.0`).
+    *   Centroid data do not retain peak width information; set `min_fwhm_ppm = 0.0` and `max_fwhm_ppm = 500.0` to avoid discarding real peaks.
+    *   Adaptive binning still works; tolerance should be based on expected mass accuracy.
+
+### Preprocessing Is Not a One-Shot Process
+The recommended parameters are a starting point, not a final answer.
+*   **Reproducibility**: Keep a record of all parameters used. The `params` dictionary can be saved as a **JSON**. This is essential for publications.
+*   **Biological Context**: If looking for rare metabolites appearing in few pixels, lower the `frequency_threshold` in binning. If interested only in high-abundance lipids, increase the SNR threshold.
+*   **Modularity**: The pipeline is modular; you can always rerun parts of it with adjusted parameters without starting from scratch.
+
+---
+*Any revision to the code, bug, error, report to julian.sierrag@icloud.com, if you we're to cite this file, remember i provide no citation and its based on internet knowledge of the web about these algorithms, so do your own research and cite properly. to cite JuliaMSI platform, please cite:*
+
+José Julián Sierra-Álvarez, Martín Orlando Camargo-Escalante, Carlos Daniel Sierra-Álvarez, Carmelo Hernández-Caricio, Juan Francisco Moreno-Luna, Isabel Buendía-Corona, Robert Winkler,
+JuliaMSI: A high-performance graphical platform for mass spectrometry imaging data analysis,
+Analytica Chimica Acta,
+Volume 1377,
+2025,
+344613,
+ISSN 0003-2670,
+https://doi.org/10.1016/j.aca.2025.344613

Project.toml

@@ -1,4 +1,5 @@
 name = "MSI_src"
+uuid = "8f395f85-3b9c-4f7f-bf7e-12345678abcd"
 authors = ["JJSA"]
 version = "0.1.0"
 
@@ -35,6 +36,7 @@ Loess = "4345ca2d-374a-55d4-8d30-97f9976e7612"
 Mmap = "a63ad114-7e13-5084-954f-fe012c677804"
 NativeFileDialog = "e1fe445b-aa65-4df4-81c1-2041507f0fd4"
 NaturalSort = "c020b1a1-e9b0-503a-9c33-f039bfc54a85"
+PackageCompiler = "9b87118b-4619-50d2-8e1e-99f35a4d4d9d"
 Parameters = "d96e819e-fc66-5662-9728-84c9c7592b0a"
 PlotlyBase = "a03496cd-edff-5a9b-9e67-9cda94a718b5"
 Printf = "de0858da-6303-5e67-8744-51eddeeeb8d7"
@@ -42,9 +44,11 @@ ProgressMeter = "92933f4c-e287-5a05-a399-4b506db050ca"
 SavitzkyGolay = "c4bf5708-b6a6-4fbe-bcd0-6850ed671584"
 Serialization = "9e88b42a-f829-5b0c-bbe9-9e923198166b"
 Setfield = "efcf1570-3423-57d1-acb7-fd33fddbac46"
+SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"
 Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"
 StatsBase = "2913bbd2-ae8a-5f71-8c99-4fb6c76f3a91"
 StipplePlotly = "ec984513-233d-481d-95b0-a3b58b97af2b"
+Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 UUIDs = "cf7118a7-6976-5b1a-9a39-7adc72f591a4"
 
 [compat]
@@ -52,9 +56,9 @@ Accessors = "0.1"
 Base64 = "1.11"
 CSV = "0.10"
 CairoMakie = "0.13"
+CodecBase = "0.3"
 ColorSchemes = "3.30"
 Colors = "0.12"
-CodecBase = "0.3"
 ContinuousWavelets = "1.1"
 DataFrames = "1.7"
 Dates = "1.11"
@@ -63,7 +67,6 @@ GLMakie = "0.11"
 Genie = "5.31"
 GenieFramework = "2.8"
 HistogramThresholding = "0.3"
-Interpolations = "0.15"
 ImageBinarization = "0.3"
 ImageComponentAnalysis = "0.2"
 ImageContrastAdjustment = "0.3"
@@ -72,6 +75,7 @@ ImageFiltering = "0.7"
 ImageMorphology = "0.4"
 ImageSegmentation = "1.9"
 Images = "0.26"
+Interpolations = "0.15"
 JLD2 = "0.6"
 JSON = "0.21"
 Libz = "1.0"
@@ -90,5 +94,12 @@ Setfield = "1.1"
 Statistics = "1.11"
 StatsBase = "0.34"
 StipplePlotly = "0.13"
+Test = "1.11.0"
 UUIDs = "1.11"
 julia = "1.11, 1.12"
+
+[extras]
+Test = "8dfed614-e60c-5e00-aba5-a33c2d43236e"
+
+[targets]
+test = ["Test"]

README.md

@@ -3,7 +3,7 @@
 https://codeberg.org/LabABI/JuliaMSI
 
 ## Local Installation
-1. Make sure you have Julia installed and at least Julia version 1.6, If not, follow the next installation guide which contains the latest files and instructions for [juliaup](https://github.com/JuliaLang/juliaup) (recommended) alternatively, you can follow the official Julia installation instructions at [julia](https://julialang.org/downloads/platform/).
+1. Make sure you have Julia installed and at least Julia version 1.6 but not higher than 1.11, If not, follow the next installation guide which contains the latest files and instructions for [juliaup](https://github.com/JuliaLang/juliaup) (recommended) alternatively, you can follow the official Julia installation instructions at [julia](https://julialang.org/downloads/platform/).
 2. Download and decompress the following repository (You can download the repository from [download](https://codeberg.org/LabABI/JuliaMSI/archive/main.zip) or locate the option manually at the top right of the page, in the box with [•••] (More operations). For a more advanced method using Git, use this link:  https://codeberg.org/LabABI/JuliaMSI.git).
 3. Ensure that you have decompressed the ZIP file and know the working directory where it is located (the folder on your computer where you saved the JuliaMSI files).
 4. Julia Version 1.11 is required. You might need to downgrade your whole Julia installation to support it.
@@ -22,7 +22,7 @@ https://codeberg.org/LabABI/JuliaMSI
    Example of a correct route:<br>
    ~/Downloads/JuliaMSI-main/juliamsi
 2. Without entering the Julia environment, launch the project in your terminal with the following command (which works for all operating systems):
-   ```
+   ```bash
    julia --threads auto --project=. start_MSI_GUI.jl
    ```
 3. After the script has finished loading, you can open a [page](http://127.0.0.1:1481/) in your browser with the web app running.
@@ -34,6 +34,47 @@ Minimum system requirements: 4 core processor, 8 GB RAM<br>
 
 JuliaMSI is a Graphical User Interface for a library of MSI tools in Julia: https://github.com/CINVESTAV-LABI/julia_mzML_imzML
 
+## Build the System Image (One-time)
+Alternatively, you can generate custom pre-compiled `.so` / `.dll` system images by running the build script in your directory. This bakes the heavy mass-spectrometry kernels and UI dependencies into a single machine-code binary for ultra-fast startup times.
+
+### 1. Build the GUI Sysimage (for full App usage)
+```bash
+julia --project=. build_sysimage.jl
+```
+This may take 5–15 minutes. The resulting `sys_msi_gui.so` (or `.dll` on Windows) file will be around **300MB–600MB** because it contains the pre-compiled machine code for your entire graphical environment.
+*Note: A sysimage built on Linux (.so) will not work on Windows. You must run the build script once on each target operating system.*
+
+Once `sys_msi_gui.so` is created in your directory, adapt your launch command:
+```bash
+julia --project=. -e 'using Pkg; Pkg.precompile()'
+# Adjust the .so extension to .dll if on Windows or .dylib if on macOS
+julia --threads auto --project=. --sysimage sys_msi_gui.so start_MSI_GUI.jl
+```
+
+### 2. Build the Headless Sysimage (for Scripts & Data Scientists)
+If you only want to run scripts or the core `MSI_src` engine without the overhead of the Genie web server and Plotly, you can build a stripped-down, high-speed system image:
+
+```bash
+julia --project=. build_sysimage.jl --headless
+```
+
+Launch your automated tests or custom processing scripts using the headless image to bypass virtually all compilation overhead:
+```bash
+julia --threads auto --project=. --sysimage sys_msi_headless.so test/test_streaming_pipeline.jl
+```
+
+### 3. Docker Option for High-Speed Processing
+For maximum portability and execution speed without clutter, we provide a `Dockerfile.headless`. This option creates a minimal, ultra-fast container that completely omits `app.jl` and the web interface. 
+
+It automatically builds and bundles the headless sysimage, giving researchers an instant-start engine ready for heavy-duty `.imzML` batch pipelines with zero JIT latency.
+
+To build and run the Docker image:
+```bash
+docker build -t juliamsi-headless -f Dockerfile.headless .
+# Run a script mapped from your local data folder:
+docker run -it -v /my/local/data:/data juliamsi-headless julia -J sys_msi_headless.so /data/my_processing_script.jl
+```
+
 ## License
 
 JuliaMSI is published under the terms of the MIT License.
@@ -41,4 +82,10 @@ JuliaMSI is published under the terms of the MIT License.
 ## Example data
 
 Robert Winkler. (2023). mzML mass spectrometry and imzML mass spectrometry imaging test data [Data set]. 
-Zenodo. <https://doi.org/10.5281/zenodo.10084132>
+Zenodo. <https://doi.org/10.5281/zenodo.10084132>
+
+## User Guides
+
+For a general guide on how to use the JuliaMSI platform, please refer to the [JuliaMSI Main Guide](JuliaMSI_MainGuide.md).
+
+For a more detailed guide on the preprocessing pipeline, please refer to the [JuliaMSI Preprocessing Guide](JuliaMSI_Preprocessing.md).

build_sysimage.jl

@@ -0,0 +1,49 @@
+# build_sysimage.jl
+using PackageCompiler
+
+headless = "--headless" in ARGS
+
+println("--- Starting Custom System Image Build ---")
+println("Mode: ", headless ? "Headless (Scripting/Docker)" : "Full GUI (Genie App)")
+println("This process can take 5-15 minutes.")
+
+ext = Sys.iswindows() ? ".dll" : Sys.isapple() ? ".dylib" : ".so"
+sysimage_name = headless ? "sys_msi_headless" * ext : "sys_msi_gui" * ext
+sysimage_path = joinpath(@__DIR__, sysimage_name)
+
+packages_to_include = Symbol[
+    :MSI_src, # <--- Bake our core library
+    :DataFrames,
+    :SparseArrays,
+    :Mmap,
+    :Libz,
+    :Serialization,
+    :Printf,
+    :JSON
+]
+
+if !headless
+    append!(packages_to_include, [
+        :Genie,
+        :GenieFramework,
+        :PlotlyBase
+    ])
+end
+
+create_sysimage(
+    packages_to_include,
+    sysimage_path = sysimage_path,
+    precompile_execution_file = "precompile_script.jl",
+    incremental = true,
+    filter_stdlibs = false
+)
+
+println("--- System Image Build Successful! ---")
+println("Created: $(sysimage_path)")
+if headless
+    println("To use in a headles script:")
+    println("julia --threads auto --project=. --sysimage $(sysimage_name) my_script.jl")
+else
+    println("To start Julia with the GUI image, use:")
+    println("julia --threads auto --project=. --sysimage $(sysimage_name) start_MSI_GUI.jl")
+end

julia_imzML_visual.jl

@@ -1,6 +1,5 @@
 # julia_imzML_visual.jl
 
-const REGISTRY_LOCK = ReentrantLock()
 
 """
     increment_image(current_image, image_list)
@@ -904,10 +903,22 @@ function warmup_init()
     end
 end
 
+"""
+    clean_registry_path(registry_path::String)
+
+Sanitizes the file path to prevent OS/file system issues, particularly on Windows,
+by standardizing separators and removing stripping whitespace.
+"""
+function clean_registry_path(registry_path::String)
+    # Standardize to forward slashes internally for uniformity, then normpath formats for OS
+    clean_path = normpath(strip(replace(registry_path, "\\" => "/")))
+    return abspath(clean_path)
+end
+
 """
     load_registry(registry_path)
 
-Loads the dataset registry from a JSON file.
+Loads the dataset registry from a JSON file, with retries for robustness on Windows.
 
 # Arguments
 - `registry_path`: Path to the `registry.json` file.
@@ -915,77 +926,114 @@ Loads the dataset registry from a JSON file.
 # Returns
 - A dictionary containing the registry data. Returns an empty dictionary if the file doesn't exist or fails to parse.
 """
-function load_registry(registry_path)
+function load_registry(registry_path::String)
     return lock(REGISTRY_LOCK) do
-        if isfile(registry_path)
+        clean_path = clean_registry_path(registry_path)
-            try
+        if isfile(clean_path)
-                JSON.parsefile(registry_path, dicttype=Dict{String,Any})
+            for attempt in 1:5
-            catch e
+                try
-                @error "Failed to parse registry.json: $e"
+                    return JSON.parsefile(clean_path, dicttype=Dict{String,Any})
-                Dict{String,Any}()
+                catch e
+                    if attempt == 5
+                        @error "Final attempt failed to parse registry.json: $(sprint(showerror, e))"
+                        return Dict{String,Any}()
+                    else
+                        @warn "Attempt $attempt to read registry.json failed, retrying in 0.2s... Error: $(sprint(showerror, e))"
+                        sleep(0.2)
+                    end
+                end
             end
-        else
-            Dict{String,Any}()
         end
+        return Dict{String,Any}()
     end
 end
 
 """
-    extract_metadata(msi_data::MSIData, source_path::String)
+    extract_metadata(loaded_data::MSIData, full_route::String)
 
-Extracts key metadata from an `MSIData` object for display.
+Extracts key metadata from an MSIData object and file path for UI display.
-
+Returns a dictionary with a "summary" key containing parameter-value pairs.
-# Arguments
-- `msi_data`: The `MSIData` object.
-- `source_path`: The path to the source data file.
-
-# Returns
-- A dictionary containing summary statistics and metadata.
 """
-function extract_metadata(msi_data::MSIData, source_path::String)
+function extract_metadata(loaded_data::MSIData, full_route::String)
-    df = msi_data.spectrum_stats_df
+    stats = []
-    if df === nothing
+    
-        # This can happen if precompute_analytics hasn't been run
+    # Basic File Info
-        # We can still return basic info
+    push!(stats, Dict("parameter" => "Filename", "value" => basename(full_route)))
-        return Dict(
+    push!(stats, Dict("parameter" => "Source Path", "value" => full_route))
-            "summary" => [
+    
-                Dict("parameter" => "File Name", "value" => basename(source_path)),
+    # Image Dimensions
-                Dict("parameter" => "Number of Spectra", "value" => length(msi_data.spectra_metadata)),
+    w, h = loaded_data.image_dims
-                Dict("parameter" => "Image Dimensions", "value" => "$(msi_data.image_dims[1]) x $(msi_data.image_dims[2])"),
+    if w > 0 && h > 0
-            ],
+        push!(stats, Dict("parameter" => "Image Dimensions", "value" => "$w x $h"))
-            "global_min_mz" => nothing,
+    else
-            "global_max_mz" => nothing
+        push!(stats, Dict("parameter" => "Image Dimensions", "value" => "N/A (Non-imaging)"))
-        )
     end
-
+    
-    summary_stats = [
+    # Spectrum Counts
-        Dict("parameter" => "File Name", "value" => basename(source_path)),
+    num_spectra = length(loaded_data.spectra_metadata)
-        Dict("parameter" => "Number of Spectra", "value" => length(msi_data.spectra_metadata)),
+    push!(stats, Dict("parameter" => "Total Spectra", "value" => string(num_spectra)))
-        Dict("parameter" => "Image Dimensions", "value" => "$(msi_data.image_dims[1]) x $(msi_data.image_dims[2])"),
+    
-        Dict("parameter" => "Global Min m/z", "value" => @sprintf("%.4f", Threads.atomic_add!(msi_data.global_min_mz, 0.0))),
+    # m/z Range
-        Dict("parameter" => "Global Max m/z", "value" => @sprintf("%.4f", Threads.atomic_add!(msi_data.global_max_mz, 0.0))),
+    min_mz, max_mz = get_global_mz_range(loaded_data)
-        Dict("parameter" => "Mean TIC", "value" => @sprintf("%.2e", mean(df.TIC))),
+    if isfinite(min_mz) && isfinite(max_mz)
-        Dict("parameter" => "Mean BPI", "value" => @sprintf("%.2e", mean(df.BPI))),
+        push!(stats, Dict("parameter" => "Global m/z Range", "value" => "$(round(min_mz, digits=4)) - $(round(max_mz, digits=4))"))
-        Dict("parameter" => "Mean # Points", "value" => @sprintf("%.1f", mean(df.NumPoints))),
+    else
-    ]
+        push!(stats, Dict("parameter" => "Global m/z Range", "value" => "Unknown"))
-
+    end
-    if hasproperty(df, :Mode)
+    
-        centroid_count = count(==(MSI_src.CENTROID), df.Mode)
+    # Instrument Metadata
-        profile_count = count(==(MSI_src.PROFILE), df.Mode)
+    if loaded_data.instrument_metadata !== nothing
-        unknown_count = count(==(MSI_src.UNKNOWN), df.Mode)
+        inst = loaded_data.instrument_metadata
-
+        push!(stats, Dict("parameter" => "Instrument Model", "value" => inst.instrument_model))
-        push!(summary_stats, Dict("parameter" => "Centroid Spectra", "value" => string(centroid_count)))
+        push!(stats, Dict("parameter" => "Polarity", "value" => titlecase(string(inst.polarity))))
-        push!(summary_stats, Dict("parameter" => "Profile Spectra", "value" => string(profile_count)))
+        push!(stats, Dict("parameter" => "Acquisition Mode", "value" => titlecase(string(inst.acquisition_mode))))
-        if unknown_count > 0
+        
-            push!(summary_stats, Dict("parameter" => "Unknown Mode Spectra", "value" => string(unknown_count)))
+        if inst.resolution !== nothing
+            push!(stats, Dict("parameter" => "Instrument Resolution", "value" => string(inst.resolution)))
         end
     end
+    
+    return Dict("summary" => stats)
+end
 
-    return Dict(
+"""
-        "summary" => summary_stats,
+    save_registry(registry_path::String, registry_data::Dict)
-        "global_min_mz" => msi_data.global_min_mz,
+
-        "global_max_mz" => msi_data.global_max_mz
+Saves the dataset registry to a JSON file atomically to prevent data loss 
-    )
+and sharing violations on Windows. Writes to a `.tmp` file first.
+"""
+function save_registry(registry_path::String, registry_data)
+    lock(REGISTRY_LOCK) do
+        clean_path = clean_registry_path(registry_path)
+        dir = dirname(clean_path)
+        if !isdir(dir)
+            mkpath(dir)
+        end
+
+        temp_path = clean_path * ".tmp"
+
+        for attempt in 1:6
+            try
+                open(temp_path, "w") do f
+                    JSON.print(f, registry_data, 4)
+                end
+                
+                # Replace final file with temp atomically. Force=true maps to ReplaceFile/MoveFileEx
+                mv(temp_path, clean_path, force=true)
+                return true
+            catch e
+                if attempt == 6
+                    @error "Final attempt failed to write to registry.json: $(sprint(showerror, e))"
+                    isfile(temp_path) && rm(temp_path, force=true)
+                    return false
+                else
+                    @warn "Attempt $attempt to write to registry.json failed (Error: $(sprint(showerror, e))), retrying in 0.2s..."
+                    sleep(0.2)
+                end
+            end
+        end
+        return false
+    end
 end
 
 """
@@ -1000,65 +1048,22 @@ Adds or updates an entry in the dataset registry JSON file.
 - `metadata`: Optional dictionary of metadata to store.
 - `is_imzML`: Boolean indicating if the source is an imzML file.
 """
-function update_registry(registry_path, dataset_name, source_path, metadata=nothing, is_imzML=false)
+function update_registry(registry_path::String, dataset_name::String, source_path::String, metadata=nothing, is_imzML=false)
     lock(REGISTRY_LOCK) do
-        registry = if isfile(registry_path)
+        registry = load_registry(registry_path)
-            try
-                JSON.parsefile(registry_path, dicttype=Dict{String,Any})
-            catch e
-                @error "Failed to parse registry.json while updating: $e"
-                Dict{String,Any}() # Start with empty if parsing fails
-            end
-        else
-            Dict{String,Any}()
-        end
 
-        # Get existing entry if it exists, otherwise create new one
         existing_entry = get(registry, dataset_name, Dict{String,Any}())
-        
-        # Start with existing data and update only the basic fields
         entry = copy(existing_entry)
         entry["source_path"] = source_path
         entry["processed_date"] = string(now())
         entry["is_imzML"] = is_imzML
         
-        # Only update metadata if provided
         if metadata !== nothing
             entry["metadata"] = metadata
         end
         
-        # Note: has_mask and mask_path are preserved from existing_entry if they exist
-        
         registry[dataset_name] = entry
-
+        save_registry(registry_path, registry)
-        try
-            open(registry_path, "w") do f
-                JSON.print(f, registry, 4)
-            end
-        catch e
-            @error "Failed to write to registry.json: $e"
-        end
-    end
-end
-
-"""
-    save_registry(registry_path, registry_data)
-
-Saves the dataset registry to a JSON file, ensuring thread-safe access.
-
-# Arguments
-- `registry_path`: Path to the `registry.json` file.
-- `registry_data`: The dictionary containing the registry data to save.
-"""
-function save_registry(registry_path, registry_data)
-    lock(REGISTRY_LOCK) do
-        try
-            open(registry_path, "w") do f
-                JSON.print(f, registry_data, 4)
-            end
-        catch e
-            @error "Failed to write to registry.json: $e"
-        end
     end
 end
 
@@ -1126,23 +1131,24 @@ function process_file_safely(file_path, masses, params, progress_message_ref, ov
         # --- Extract metadata ---
         metadata = extract_metadata(local_msi_data, file_path)
 
-        # --- Save Slices ---
+        # --- Save Slices (Parallelized) ---
         mkpath(output_dir)
-
+        
-        for (mass_idx, mass) in enumerate(masses)
+        progress_lock = ReentrantLock()
-            progress_message_ref = "File $(params.fileIdx)/$(params.nFiles): Saving slice for m/z=$mass"
+        
-
+        # Parallelize over the requested masses to fully utilize multi-core CPUs
+        Threads.@threads for mass in masses
             slice = slice_dict[mass]
             text_nmass = replace(string(mass), "." => "_")
             bitmap_filename = params.triqE ? "TrIQ_$(text_nmass).bmp" : "MSI_$(text_nmass).bmp"
             colorbar_filename = params.triqE ? "colorbar_TrIQ_$(text_nmass).png" : "colorbar_MSI_$(text_nmass).png"
 
+            local sliceQuant, bounds
             if all(iszero, slice)
                 sliceQuant = zeros(UInt8, size(slice))
                 bounds = (0.0, 1.0)  # Default bounds for empty slice
-                @warn "No intensity data for m/z = $mass in $(dataset_name)"
             else
-                # Get both quantized data AND bounds in one call
+                # TrIQ and quantization are computationally intensive and now run in parallel
                 if params.triqE
                     sliceQuant, bounds = TrIQ(slice, params.colorL, params.triqP, mask_matrix=mask_matrix_for_triq)
                 else
@@ -1151,18 +1157,24 @@ function process_file_safely(file_path, masses, params, progress_message_ref, ov
                 
                 if params.medianF
                     sliceQuant = round.(UInt8, median_filter(sliceQuant))
-                    # Note: bounds remain the same after median filter
                 end
             end
 
+            # Disk I/O (Saving BMP)
             save_bitmap(joinpath(output_dir, bitmap_filename), sliceQuant, ViridisPalette)
+            
             if !all(iszero, slice)
-                # Now pass the precomputed bounds to colorbar generation
+                # CairoMakie colorbar generation is now also parallelized
                 generate_colorbar_image(slice, params.colorL, joinpath(output_dir, colorbar_filename), 
                                       bounds; use_triq=params.triqE, triq_prob=params.triqP, mask_path=mask_path)
             end
+            
+            lock(progress_lock) do
+                progress_message_ref = "File $(params.fileIdx)/$(params.nFiles): Saved slice for m/z=$mass"
+            end
         end
 
+
         is_imzML = local_msi_data.source isa ImzMLSource
         update_registry(params.registry, dataset_name, file_path, metadata, is_imzML)
         return (true, "")
@@ -1182,24 +1194,14 @@ function process_file_safely(file_path, masses, params, progress_message_ref, ov
     end
 end
 
-function update_registry_mask_fields(registry_path, dataset_name, has_mask, mask_path)
+function update_registry_mask_fields(registry_path::String, dataset_name::String, has_mask::Bool, mask_path::String)
     lock(REGISTRY_LOCK) do
-        registry = if isfile(registry_path)
+        registry = load_registry(registry_path)
-            try
-                JSON.parsefile(registry_path, dicttype=Dict{String,Any})
-            catch e
-                @error "Failed to parse registry.json while updating mask fields: $e"
-                Dict{String,Any}()
-            end
-        else
-            Dict{String,Any}()
-        end
         
         if haskey(registry, dataset_name)
             registry[dataset_name]["has_mask"] = has_mask
             registry[dataset_name]["mask_path"] = mask_path
         else
-            # Create new entry if it doesn't exist
             registry[dataset_name] = Dict{String,Any}(
                 "has_mask" => has_mask,
                 "mask_path" => mask_path,
@@ -1209,13 +1211,7 @@ function update_registry_mask_fields(registry_path, dataset_name, has_mask, mask
             )
         end
 
-        try
+        save_registry(registry_path, registry)
-            open(registry_path, "w") do f
-                JSON.print(f, registry, 4)
-            end
-        catch e
-            @error "Failed to write to registry.json: $e"
-        end
     end
 end
 

mask.jl

@@ -9,10 +9,12 @@ using Statistics, NaturalSort, LinearAlgebra, StipplePlotly
 using Base.Filesystem: mv
 
 using MSI_src
-using .MSI_src: MSIData, process_image_pipeline
+using .MSI_src: MSIData, process_image_pipeline, REGISTRY_LOCK
 
 # Plot Handling
-include("./julia_imzML_visual.jl")
+if !@isdefined(increment_image)
+    include("./julia_imzML_visual.jl")
+end
 
 # Image Processing Pipeline
 using ImageBinarization
@@ -799,9 +801,7 @@ end
                 end
                 
                 # Save the updated registry
-                open(reg_path, "w") do f
+                save_registry(reg_path, registry)
-                    JSON.print(f, registry, 4)
-                end
                 
                 @info "Registry updated with mask: $(final_mask_name)"
                 
@@ -874,9 +874,7 @@ end
 
                 if !isempty(new_folders) || !isempty(removed_folders)
                     println("Registry changed, saving...")
-                    open(reg_path, "w") do f
+                    save_registry(reg_path, registry)
-                        JSON.print(f, registry, 4)
-                    end
                 end
                 
                 all_folders = sort(collect(keys(registry)), lt=natural)

precompile_script.jl

@@ -0,0 +1,173 @@
+# precompile_script.jl
+# This script exercises the core MSI processing kernels (imzML and mzML) 
+# to ensure they are precompiled into the custom system image.
+
+using MSI_src
+using Mmap
+using Base64
+using PlotlyBase
+
+# --- DYNAMIC WARM-UP WORKLOAD --- #
+
+function create_warmup_datasets(dir, T_mz::Type, T_int::Type)
+    # 1. Create minimal imzML/ibd
+    imzml_path = joinpath(dir, "warmup_$(T_mz)_$(T_int).imzML")
+    ibd_path = joinpath(dir, "warmup_$(T_mz)_$(T_int).ibd")
+    
+    n_points = 10
+    n_pixels = 1
+    mz_bytes = n_points * sizeof(T_mz)
+    int_bytes = n_points * sizeof(T_int)
+    
+    # Write some dummy binary data (m/z must be sorted for validation)
+    open(ibd_path, "w") do f
+        write(f, sort(rand(T_mz, n_points)))
+        write(f, rand(T_int, n_points))
+    end
+    
+    # Minimal but VALID imzML structure that passes axes_config_img
+    imzml_content = """<?xml version="1.0" encoding="utf-8"?>
+<mzML xmlns="http://psi.hupo.org/ms/mzML" version="1.1.0">
+  <cvList count="2">
+    <cv id="MS" fullName="Proteomics Standards Initiative Mass Spectrometry Ontology" version="3.30.0" URI="http://psidev.cvs.sourceforge.net/*checkout*/psidev/ms/mzML/ontology/psi-ms.obo"/>
+    <cv id="IMS" fullName="Imaging MS Ontology" version="0.9.1" URI="http://www.imzml.org/ontology/imzML1.1.0.obo"/>
+  </cvList>
+  <fileDescription>
+    <fileContent>
+      <cvParam cvRef="IMS" accession="IMS:1000031" name="processed" value=""/>
+      <cvParam cvRef="IMS" accession="IMS:1000080" name="universally unique identifier" value="00000000-0000-0000-0000-000000000000"/>
+      <cvParam cvRef="IMS" accession="IMS:1000091" name="ibd MD5 HTTP" value="00000000000000000000000000000000"/>
+    </fileContent>
+  </fileDescription>
+  <referenceableParamGroupList count="2">
+    <referenceableParamGroup id="mz_array_settings">
+      <cvParam cvRef="MS" accession="MS:1000514" name="m/z array" value=""/>
+      <cvParam cvRef="MS" accession="MS:$(T_mz == Float64 ? "1000523" : "1000521")" name="$(T_mz == Float64 ? "64-bit float" : "32-bit float")" value=""/>
+      <cvParam cvRef="MS" accession="MS:1000574" name="zlib compression" value=""/>
+    </referenceableParamGroup>
+    <referenceableParamGroup id="intensity_array_settings">
+      <cvParam cvRef="MS" accession="MS:1000515" name="intensity array" value=""/>
+      <cvParam cvRef="MS" accession="MS:$(T_int == Float64 ? "1000523" : "1000521")" name="$(T_int == Float64 ? "64-bit float" : "32-bit float")" value=""/>
+      <cvParam cvRef="MS" accession="MS:1000574" name="zlib compression" value=""/>
+    </referenceableParamGroup>
+  </referenceableParamGroupList>
+  <run id="run_1">
+    <spectrumList count="1">
+      <spectrum index="0" id="scan=1" defaultArrayLength="$n_points">
+        <cvParam cvRef="IMS" accession="IMS:1000050" name="coordinate x" value="1"/>
+        <cvParam cvRef="IMS" accession="IMS:1000051" name="coordinate y" value="1"/>
+        <binaryDataArrayList count="2">
+          <binaryDataArray encodedLength="0">
+            <referenceableParamGroupRef ref="mz_array_settings"/>
+            <cvParam cvRef="IMS" accession="IMS:1000102" name="external data" value=""/>
+            <cvParam cvRef="IMS" accession="IMS:1000103" name="external offset" value="0"/>
+            <cvParam cvRef="IMS" accession="IMS:1000104" name="external array length" value="$n_points"/>
+            <cvParam cvRef="IMS" accession="IMS:1000101" name="external encoded length" value="$mz_bytes"/>
+            <binary/>
+          </binaryDataArray>
+          <binaryDataArray encodedLength="0">
+            <referenceableParamGroupRef ref="intensity_array_settings"/>
+            <cvParam cvRef="IMS" accession="IMS:1000102" name="external data" value=""/>
+            <cvParam cvRef="IMS" accession="IMS:1000103" name="external offset" value="$mz_bytes"/>
+            <cvParam cvRef="IMS" accession="IMS:1000104" name="external array length" value="$n_points"/>
+            <cvParam cvRef="IMS" accession="IMS:1000101" name="external encoded length" value="$int_bytes"/>
+            <binary/>
+          </binaryDataArray>
+        </binaryDataArrayList>
+      </spectrum>
+    </spectrumList>
+  </run>
+</mzML>
+"""
+    write(imzml_path, imzml_content)
+
+    # 2. Create minimal mzML (Base64)
+    mzml_path = joinpath(dir, "warmup_$(T_mz)_$(T_int).mzML")
+    b64_mz = base64encode(sort(rand(T_mz, n_points)))
+    b64_int = base64encode(rand(T_int, n_points))
+    
+    mzml_content = """<?xml version="1.0" encoding="utf-8"?>
+<mzML xmlns="http://psi.hupo.org/ms/mzML" version="1.1.0">
+  <run id="run_1">
+    <spectrumList count="1">
+      <spectrum index="0" id="scan=1" defaultArrayLength="$n_points">
+        <binaryDataArrayList count="2">
+          <binaryDataArray encodedLength="$(length(b64_mz))">
+            <cvParam cvRef="MS" accession="MS:1000514" name="m/z array" value=""/>
+            <cvParam cvRef="MS" accession="MS:$(T_mz == Float64 ? "1000523" : "1000521")" name="" value=""/>
+            <binary>$b64_mz</binary>
+          </binaryDataArray>
+          <binaryDataArray encodedLength="$(length(b64_int))">
+            <cvParam cvRef="MS" accession="MS:1000515" name="intensity array" value=""/>
+            <cvParam cvRef="MS" accession="MS:$(T_int == Float64 ? "1000523" : "1000521")" name="" value=""/>
+            <binary>$b64_int</binary>
+          </binaryDataArray>
+        </binaryDataArrayList>
+      </spectrum>
+    </spectrumList>
+  </run>
+  <indexList count="1"><index name="spectrum"><offset idRef="scan=1">PLACEHOLDER</offset></index></indexList>
+  <indexListOffset>ID_OFFSET</indexListOffset>
+</mzML>
+"""
+    spec_start = findfirst("<spectrum", mzml_content).start - 1
+    index_start = findfirst("<indexList", mzml_content).start - 1
+    mzml_content = replace(mzml_content, "PLACEHOLDER" => string(spec_start))
+    mzml_content = replace(mzml_content, "ID_OFFSET" => string(index_start))
+    
+    write(mzml_path, mzml_content)
+    
+    return imzml_path, mzml_path
+end
+
+println("Starting robust precompilation warmup (imzML + mzML)...")
+
+try
+    mktempdir() do tmp_dir
+        for (T_mz, T_int) in [(Float32, Float32), (Float64, Float32)]
+            println("Exercising pathways for $T_mz and $T_int...")
+            imzml_path, mzml_path = create_warmup_datasets(tmp_dir, T_mz, T_int)
+            
+            try
+                # Exercise imzML pathway
+                imzml_data = MSI_src.load_imzml_lazy(imzml_path; use_mmap=true)
+                MSI_src.get_multiple_mz_slices(imzml_data, [10.0, 20.0], 0.1)
+                
+                # Exercise Sprint 2 Streaming Pipeline
+                config = MSI_src.PipelineConfig(
+                    steps=[
+                        MSI_src.StreamingStep(:baseline_correction, Dict(:method => :snip, :iterations => 5)),
+                        MSI_src.StreamingStep(:normalization, Dict(:method => :tic)),
+                        MSI_src.StreamingStep(:peak_picking, Dict(:method => :profile, :snr_threshold => 3.0))
+                    ],
+                    num_bins=2000
+                )
+                try
+                    MSI_src.process_dataset!(imzml_data, config)
+                catch
+                    # Ignore matrix mismatch errors from tiny random data 
+                end
+                
+                # Exercise mzML pathway
+                mzml_data = MSI_src.load_mzml_lazy(mzml_path)
+                MSI_src.GetSpectrum(mzml_data, 1)
+            catch e
+                @warn "Warmup failed for $T_mz/$T_int: $e"
+            end
+        end
+
+        # --- COMMON KERNELS ---
+        println("Precompiling image processing kernels...")
+        dummy_slice = rand(Float32, 4, 4)
+        MSI_src.TrIQ(dummy_slice, 256, 1.0)
+        MSI_src.save_bitmap(joinpath(tmp_dir, "d.bmp"), zeros(UInt8, 4, 4), MSI_src.ViridisPalette)
+        
+        # Plotly Warmup
+        p = PlotlyBase.Plot(PlotlyBase.scatter(x=1:2, y=[1,2]))
+    end
+catch e
+    if !(e isa InterruptException)
+        @warn "Global Warmup failed: $e"
+    end
+end
+println("Precompilation workload completed.")

src/BloomFilters.jl

@@ -20,21 +20,7 @@ mutable struct BloomFilter{T}
     count::Int
 end
 
-"""
-    StreamingBloomFilter
 
-A memory-efficient Bloom filter that doesn't store all values at once.
-"""
-mutable struct StreamingBloomFilter
-    bits::BitVector
-    size::Int
-    hash_count::Int
-    seed::UInt64
-    count::Int
-    # Streaming state
-    current_chunk::Vector{Int}
-    chunk_size::Int
-end
 
 """
     BloomFilter{T}(expected_elements::Int, false_positive_rate::Float64=0.01; kwargs...) -> Return type
@@ -51,6 +37,9 @@ Creates a Bloom filter optimized for the expected number of elements and desired
 
 - `seed::Union{UInt32,UInt64}`: Keyword description
     (**Default**: `0x12345678`)
+
+# Returns
+- A `BloomFilter{T}`.
 """
 function BloomFilter{T}(expected_elements::Int, false_positive_rate::Float64=0.01; seed::Union{UInt32,UInt64}=0x12345678) where T
     # Convert seed to UInt64 for consistency
@@ -73,6 +62,9 @@ Calculates the optimal number of bits for a Bloom filter
 
 - `n::Int`: Expected number of elements
 - `p::Float64`: Desired false positive rate
+
+# Returns
+- An `Int` to determine the optimal number of bits for a Bloom filter.
 """
 function optimal_bit_size(n::Int, p::Float64)::Int
     if p <= 0.0 || p >= 1.0
@@ -91,7 +83,10 @@ Calculates the optimal number of hash functions for a Bloom filter.
 # Arguments
 
 - `n::Int`: Expected number of elements
-- `p::Float64`: Desired false positive rate
+- `m::Int`: Number of bits in the filter
+
+# Returns
+- An `Int` to determine the optimal number of hash functions for a Bloom filter.
 """
 function optimal_hash_count(n::Int, m::Int)::Int
     if n <= 0 || m <= 0
@@ -106,6 +101,16 @@ end
     hash_functions(item::T, count::Int, size::Int, seed::UInt64) -> Vector{Int}
 
 Generates multiple hash values for an item using double hashing technique.
+
+# Arguments
+
+- `item::T`: Argument description
+- `count::Int`: Argument description
+- `size::Int`: Argument description
+- `seed::UInt64`: Argument description
+
+# Returns
+- A `Vector{Int}` of hash values.
 """
 function hash_functions(item::T, count::Int, size::Int, seed::UInt64) where T
     # Use Julia's built-in hash with different seeds
@@ -129,7 +134,7 @@ end
 
 Adds an element to the Bloom filter.
 
-# Arguments
+# Arguments:
 
 - `bf::BloomFilter{T}`: Argument description
 - `item::T`: Argument description
@@ -150,6 +155,11 @@ end
 
 Checks whether an element is possibly in the Bloom filter.
 Returns `true` if the element might be in the set, `false` if it's definitely not.
+
+# Arguments:
+
+- `item::T`: Argument description
+- `bf::BloomFilter{T}`: Argument description
 """
 function Base.in(item::T, bf::BloomFilter{T})::Bool where T
     hashes = hash_functions(item, bf.hash_count, bf.size, bf.seed)
@@ -166,6 +176,11 @@ end
     contains(bf::BloomFilter{T}, item::T) -> Bool
 
 Alias for `in()` for compatibility with your existing code.
+
+# Arguments:
+
+- `bf::BloomFilter{T}`: Argument description
+- `item::T`: Argument description
 """
 contains(bf::BloomFilter{T}, item::T) where T = item in bf
 
@@ -173,6 +188,11 @@ contains(bf::BloomFilter{T}, item::T) where T = item in bf
     add!(bf::BloomFilter{T}, item::T)
 
 Alias for `push!()` for compatibility with your existing code.
+
+# Arguments:
+
+- `bf::BloomFilter{T}`: Argument description
+- `item::T`: Argument description
 """
 add!(bf::BloomFilter{T}, item::T) where T = push!(bf, item)
 
@@ -180,6 +200,10 @@ add!(bf::BloomFilter{T}, item::T) where T = push!(bf, item)
     false_positive_rate(bf::BloomFilter) -> Float64
 
 Estimates the current false positive rate of the Bloom filter.
+
+# Arguments:
+
+- `bf::BloomFilter`: Argument description
 """
 function false_positive_rate(bf::BloomFilter)::Float64
     if bf.count == 0

src/Common.jl

@@ -1,111 +1,113 @@
-# src/Common.jl - Updated with BloomFilter
 using Base.Threads
+using Mmap
+
+# POSIX madvise constants
+const MADV_NORMAL = 0
+const MADV_RANDOM = 1
+const MADV_SEQUENTIAL = 2
+const MADV_WILLNEED = 3
+const MADV_DONTNEED = 4
+
+"""
+    posix_madvise(buffer::AbstractArray, advice::Integer)
+
+A safe wrapper for the OS `madvise` system call. Signals the kernel about the 
+access pattern for a memory-mapped region. Currently supports Linux/Unix systems.
+"""
+function posix_madvise(buffer::AbstractArray, advice::Integer)
+    @static if Sys.isunix()
+        try
+            ptr = pointer(buffer)
+            len = sizeof(buffer)
+            # ccall(:madvise, return_type, (arg_types...), args...)
+            ret = ccall(:madvise, Int32, (Ptr{Cvoid}, Csize_t, Int32), ptr, len, Int32(advice))
+            return ret == 0
+        catch
+            return false
+        end
+    else
+        return false # Not supported on this OS
+    end
+end
 
 # --- Buffer Pooling ---
-"""
-    BufferPool
 
-A thread-safe pool of `Vector{UInt8}` buffers, categorized by their size.
-This helps reduce memory allocations and garbage collection overhead by reusing buffers.
-
-# Fields
-- `lock`: A `ReentrantLock` to ensure thread-safe access to the pool.
-- `buffers`: A dictionary mapping buffer size (in bytes) to a list of available buffers of that size.
-"""
-mutable struct BufferPool
-    lock::ReentrantLock
-    buffers::Dict{Int, Vector{Vector{UInt8}}}
-
-    function BufferPool()
-        new(ReentrantLock(), Dict{Int, Vector{Vector{UInt8}}}())
-    end
-end
-
-"""
-    get_buffer!(pool::BufferPool, size::Int) -> Vector{UInt8}
-
-Retrieves a `Vector{UInt8}` buffer of at least `size` bytes from the pool.
-If no suitable buffer is available, a new one is allocated.
-
-# Arguments
-- `pool`: The `BufferPool` to retrieve the buffer from.
-- `size`: The minimum desired size of the buffer in bytes.
-
-# Returns
-- A `Vector{UInt8}` buffer.
-"""
-function get_buffer!(pool::BufferPool, size::Int)::Vector{UInt8}
-    lock(pool.lock) do
-        # Find an existing buffer that is large enough
-        for (buffer_size, buffer_list) in pool.buffers
-            if buffer_size >= size && !isempty(buffer_list)
-                buffer = pop!(buffer_list)
-                # Resize if necessary (and if it's significantly larger than needed)
-                if length(buffer) > size * 2 # Heuristic: if buffer is more than twice the requested size
-                    return Vector{UInt8}(undef, size) # Allocate new smaller buffer
-                else
-                    return buffer
-                end
-            end
-        end
-        # No suitable buffer found, allocate a new one
-        return Vector{UInt8}(undef, size)
-    end
-end
-
-"""
-    release_buffer!(pool::BufferPool, buffer::Vector{UInt8})
-
-Returns a `Vector{UInt8}` buffer to the pool for future reuse.
-
-# Arguments
-- `pool`: The `BufferPool` to return the buffer to.
-- `buffer`: The `Vector{UInt8}` buffer to release.
-"""
-function release_buffer!(pool::BufferPool, buffer::Vector{UInt8})
-    lock(pool.lock) do
-        buffer_size = length(buffer)
-        if !haskey(pool.buffers, buffer_size)
-            pool.buffers[buffer_size] = Vector{Vector{UInt8}}()
-        end
-        push!(pool.buffers[buffer_size], buffer)
-    end
-    return nothing
-end
 
 
 """
     SimpleBufferPool
 
 A dramatically simplified buffer pool that avoids complex locking.
+
+# Arguments:
+- `buffers::Dict{Int, Vector{Vector{UInt8}}}`: A dictionary mapping buffer size (in bytes) to a list of available buffers of that size.
+- `max_pool_size::Int`: The maximum number of buffers to keep in the pool.
+- `lock::ReentrantLock`: A lock to ensure thread-safe access to the pool.
 """
 mutable struct SimpleBufferPool
     buffers::Dict{Int, Vector{Vector{UInt8}}}
     max_pool_size::Int
+    lock::ReentrantLock  # Add lock for thread safety
 end
 
-SimpleBufferPool() = SimpleBufferPool(Dict{Int, Vector{Vector{UInt8}}}(), 50)
+"""
+    SimpleBufferPool()
 
+Creates a new `SimpleBufferPool` with a default maximum pool size of 50.
+"""
+SimpleBufferPool() = SimpleBufferPool(Dict{Int, Vector{Vector{UInt8}}}(), 50, ReentrantLock())
+
+"""
+    get_buffer!(pool::SimpleBufferPool, size::Int) -> Vector{UInt8}
+
+Retrieves a `Vector{UInt8}` buffer of at least `size` bytes from the pool.
+If no suitable buffer is available, a new one is allocated.
+
+# Arguments:
+- `pool`: The `SimpleBufferPool` to retrieve the buffer from.
+- `size`: The minimum desired size of the buffer in bytes.
+
+# Returns:
+- A `Vector{UInt8}` buffer.
+"""
 function get_buffer!(pool::SimpleBufferPool, size::Int)::Vector{UInt8}
-    # Check for existing buffers of exact size first
+    buf = lock(pool.lock) do
-    if haskey(pool.buffers, size) && !isempty(pool.buffers[size])
+        # Check for existing buffers of exact size first
-        return pop!(pool.buffers[size])
+        if haskey(pool.buffers, size) && !isempty(pool.buffers[size])
+            return pop!(pool.buffers[size])
+        end
+        return nothing
     end
     
-    # No suitable buffer found, allocate new one
+    if buf !== nothing
+        return buf
+    end
+    
+    # No suitable buffer found, allocate new one (outside lock to reduce contention)
     return Vector{UInt8}(undef, size)
 end
 
+"""
+    release_buffer!(pool::SimpleBufferPool, buffer::Vector{UInt8})
+
+Returns a `Vector{UInt8}` buffer to the pool for future reuse.
+
+# Arguments:
+- `pool`: The `SimpleBufferPool` to return the buffer to.
+- `buffer`: The `Vector{UInt8}` buffer to release.
+"""
 function release_buffer!(pool::SimpleBufferPool, buffer::Vector{UInt8})
     size = length(buffer)
     
-    if !haskey(pool.buffers, size)
+    lock(pool.lock) do
-        pool.buffers[size] = Vector{Vector{UInt8}}()
+        if !haskey(pool.buffers, size)
-    end
+            pool.buffers[size] = Vector{Vector{UInt8}}()
-    
+        end
-    # Limit pool size to prevent memory bloat
+        
-    if length(pool.buffers[size]) < pool.max_pool_size
+        # Limit pool size to prevent memory bloat
-        push!(pool.buffers[size], buffer)
+        if length(pool.buffers[size]) < pool.max_pool_size
+            push!(pool.buffers[size], buffer)
+        end
     end
     # If pool is full, let buffer get GC'd
 end
@@ -127,15 +129,26 @@ struct InvalidSpectrumError <: MSIError
 end
 
 # --- Shared Constants ---
-const DEFAULT_CACHE_SIZE = 100
+const DEFAULT_CACHE_SIZE = 100 # Default cache size for spectra
-const DEFAULT_NUM_BINS = 2000
+const DEFAULT_NUM_BINS = 2000 # Default number of bins for spectra
 const DEFAULT_BLOOM_FILTER_SIZE = 10000  # Default size for empty spectra
-const DEFAULT_FALSE_POSITIVE_RATE = 0.01
+const DEFAULT_FALSE_POSITIVE_RATE = 0.01 # Default false positive rate for bloom filters
 
 """
     validate_spectrum_data(mz, intensity, id)
 
 Checks a spectrum for common data integrity issues.
+
+# Arguments:
+- `mz`: The m/z values of the spectrum.
+- `intensity`: The intensity values of the spectrum.
+- `id`: The ID of the spectrum.
+
+# Returns:
+- `true` if the spectrum is valid.
+
+# Throws:
+- `InvalidSpectrumError` if the spectrum is invalid.
 """
 function validate_spectrum_data(mz::AbstractVector, intensity::AbstractVector, id)
     if length(mz) != length(intensity)
@@ -158,6 +171,12 @@ end
 
 Memory-optimized version that processes data in chunks to reduce temporary allocations
 and avoid holding large intermediate arrays.
+
+# Arguments:
+- `mz`: The m/z values of the spectrum.
+
+# Returns:
+- A `BloomFilter{Int}` for the spectrum.
 """
 function create_bloom_filter_for_spectrum(mz::AbstractVector{<:Real})::BloomFilter{Int}
     if isempty(mz)
@@ -188,6 +207,9 @@ end
     empty_bloom_filter() -> BloomFilter{Float64}
 
 Creates an empty Bloom filter for spectra with no data.
+
+# Returns:
+- An empty `BloomFilter{Float64}`.
 """
 function empty_bloom_filter()::BloomFilter{Int}
     return BloomFilter{Int}(size=DEFAULT_BLOOM_FILTER_SIZE, hash_count=3)
@@ -197,6 +219,12 @@ end
     AtomicFlag
 
 A thread-safe boolean flag using atomic operations.
+
+# Fields:
+- `value::Base.Threads.Atomic{Int}`: The atomic value of the flag.
+
+# Arguments:
+- `value::Base.Threads.Atomic{Int}`: The atomic value of the flag.
 """
 struct AtomicFlag
     value::Base.Threads.Atomic{Int}

src/MSI_src.jl

@@ -19,7 +19,22 @@ export OpenMSIData,
        get_global_mz_range,
        MSIData,
        _iterate_spectra_fast,
-       validate_spectrum
+       validate_spectrum,
+       get_mz_slice,
+       REGISTRY_LOCK,
+       SpectrumMetadata,
+       SpectrumAsset,
+       SpectrumMode,
+       CENTROID,
+       PROFILE,
+       MzMLSource,
+       ImzMLSource,
+       SimpleBufferPool,
+       get_buffer!,
+       release_buffer!
+
+# Define shared registry lock
+const REGISTRY_LOCK = ReentrantLock()
 
 # Export the public preprocessing & precalculations API
 export run_preprocessing_analysis,
@@ -56,9 +71,21 @@ export apply_baseline_correction,
        apply_intensity_transformation,
        save_feature_matrix
 
+# Sprint 2: Streaming Pipeline API
+export process_dataset!,
+       PipelineConfig,
+       StreamingStep,
+       normalize_inplace!,
+       transform_inplace!,
+       smooth_inplace!,
+       baseline_subtract_inplace!,
+       detect_peaks_streaming,
+       calibrate_inplace!
+
 # Include all source files directly into the main module
 include("BloomFilters.jl")
 include("Common.jl")
+include("ResourcePool.jl")
 include("MSIData.jl")
 include("ParserHelpers.jl")
 include("mzML.jl")
@@ -68,6 +95,8 @@ include("Preprocessing.jl")
 include("ImageProcessing.jl")
 include("Precalculations.jl")
 include("PreprocessingPipeline.jl")
+include("StreamingKernels.jl")
+include("StreamingPipeline.jl")
 
 using Setfield # For immutable struct updates
 
@@ -80,15 +109,28 @@ using Setfield # For immutable struct updates
 Opens a .mzML or .imzML file and prepares it for data access.
 
 This is the main entry point for the new data access API.
+
+# Arguments
+
+- `filepath::String`: Path to the .mzML or .imzML file to open.
+- `cache_size::Int`: Number of spectra to cache in memory.
+- `spectrum_type_map::Union{Dict{Int, Symbol}, Nothing}`: Optional mapping of spectrum indices to types.
+
+# Returns
+- An `MSIData` object.
 """
 function OpenMSIData(filepath::String; cache_size=300, spectrum_type_map::Union{Dict{Int, Symbol}, Nothing}=nothing)
+    # Apply standard path normalization for cross-platform compatibility
+    # Ensure Windows backslashes are converted to OS-native separators
+    norm_filepath = normpath(replace(filepath, "\\" => "/"))
+    
     local msi_data
-    if endswith(lowercase(filepath), ".mzml")
+    if endswith(lowercase(norm_filepath), ".mzml")
-        msi_data = load_mzml_lazy(filepath, cache_size=cache_size)
+        msi_data = load_mzml_lazy(norm_filepath, cache_size=cache_size)
-    elseif endswith(lowercase(filepath), ".imzml")
+    elseif endswith(lowercase(norm_filepath), ".imzml")
-        msi_data = load_imzml_lazy(filepath, cache_size=cache_size)
+        msi_data = load_imzml_lazy(norm_filepath, cache_size=cache_size)
     else
-        error("Unsupported file type: $filepath. Please provide a .mzML or .imzML file.")
+        error("Unsupported file type: $norm_filepath. Please provide a .mzML or .imzML file.")
     end
 
     # Apply spectrum type map if provided

src/MzmlConverter.jl

@@ -842,6 +842,16 @@ Main workflow function to convert a .mzML file to an .imzML file.
 * `target_file`: Path for the output .imzML file (the .ibd will be named accordingly).
 * `img_width`: width dimention for the creation of the x axis
 * `img_height`: height dimention for the creation of the y axis
+
+# Returns
+- `true` if the conversion was successful.
+- `false` if the conversion failed.
+
+# Example
+
+```julia
+ImportMzmlFile("test.mzML", "test.txt", "test.imzML")
+```
 """
 function ImportMzmlFile(source_file::String, sync_file::String, target_file::String; img_width::Int=0, img_height::Int=0)
     if !isfile(source_file)

src/ParserHelpers.jl

@@ -28,8 +28,18 @@ end
 
 Reads a stream line-by-line until a line matches the provided regex.
 
+# Arguments
+- `stream::IO`: The stream to parse the metadata from.
+- `regex::Regex`: The regex to match.
+
 # Returns
 - A `RegexMatch` object if a match is found, otherwise throws an error.
+
+# Example
+
+```julia
+find_tag(open("test.mzML"), r"<spectrum")
+```
 """
 function find_tag(stream, regex::Regex)
     while !eof(stream)
@@ -47,8 +57,18 @@ end
 
 Retrieves an attribute's value from an XML tag string.
 
+# Arguments
+- `source::String`: The string to parse the attribute from.
+- `tag::String`: The tag to match.
+
 # Returns
 - A `RegexMatch` object containing the attribute and its value.
+
+# Example
+
+```julia
+get_attribute("<spectrum index=\"1\">", "index")
+```
 """
 function get_attribute(source::AbstractString, tag::String = "([^=]+)")
     # Construct the regex pattern string
@@ -65,10 +85,16 @@ struct. It reads accessions to determine the data format (e.g., `Float32`),
 compression status (`zlib`), and axis type (m/z vs. intensity).
 
 # Arguments
-- `stream`: An IO stream positioned at the start of the `cvParam` block.
+- `stream::IO`: An IO stream positioned at the start of the `cvParam` block.
 
 # Returns
 - A `SpecDim` struct populated with the parsed configuration.
+
+# Example
+
+```julia
+configure_spec_dim(open("test.mzML"))
+```
 """
 function configure_spec_dim(stream)
     axis = SpecDim(Float64, false, 1, 0, UNKNOWN)  # Add UNKNOWN as default mode

src/Precalculations.jl

@@ -1,4 +1,5 @@
 using StatsBase # For mean, std, median, quantile, mad
+using Base.Threads # For atomic counters and locking
 
 # =============================================================================
 # 7) Spatial & Advanced Processing (Stubs & New Functions)
@@ -80,11 +81,12 @@ function analyze_mass_accuracy(
     println("Analyzing mass accuracy for $(length(spectrum_indices)) spectra...")
 
     all_ppm_errors = Float64[]
-    total_matched_peaks = 0
+    total_matched_peaks = Base.Threads.Atomic{Int}(0)
-    total_spectra_processed = 0
+    total_spectra_processed = Base.Threads.Atomic{Int}(0)
+    results_lock = ReentrantLock()
 
     _iterate_spectra_fast(msi_data, spectrum_indices) do idx, mz, intensity
-        total_spectra_processed += 1
+        Base.Threads.atomic_add!(total_spectra_processed, 1)
         if !validate_spectrum(mz, intensity)
             @warn "Spectrum $idx is invalid, skipping mass accuracy analysis for it."
             return
@@ -106,8 +108,10 @@ function analyze_mass_accuracy(
             end
 
             if best_matched_peak_mz !== nothing
-                push!(all_ppm_errors, min_ppm_error)
+                lock(results_lock) do
-                total_matched_peaks += 1
+                    push!(all_ppm_errors, min_ppm_error)
+                end
+                Base.Threads.atomic_add!(total_matched_peaks, 1)
             end
         end
     end
@@ -140,8 +144,8 @@ function analyze_mass_accuracy(
         std_ppm_error = std_err,
         min_ppm_error = min_err,
         max_ppm_error = max_err,
-        total_matched_peaks = total_matched_peaks,
+        total_matched_peaks = total_matched_peaks[],
-        total_spectra_analyzed = total_spectra_processed,
+        total_spectra_analyzed = total_spectra_processed[],
         ppm_error_distribution = all_ppm_errors
     )
 end
@@ -318,21 +322,29 @@ function analyze_instrument_characteristics(msi_data::MSIData; sample_indices::A
         return results
     end
     
+    results_lock = ReentrantLock()
+    
     _iterate_spectra_fast(msi_data, sample_indices) do idx, mz, intensity
         # Record spectrum mode
-        push!(spectrum_modes, msi_data.spectra_metadata[idx].mode)
+        lock(results_lock) do
+            push!(spectrum_modes, msi_data.spectra_metadata[idx].mode)
+        end
         
         # Calculate m/z step statistics (for profile data)
         if length(mz) > 1 && msi_data.spectra_metadata[idx].mode == PROFILE
             steps = diff(mz)
             if !isempty(steps)
-                push!(mz_step_sizes, mean(steps))
+                lock(results_lock) do
+                    push!(mz_step_sizes, mean(steps))
+                end
             end
         end
         
         # Record intensity range
         if !isempty(intensity)
-            push!(intensity_ranges, (minimum(intensity), maximum(intensity)))
+            lock(results_lock) do
+                push!(intensity_ranges, (minimum(intensity), maximum(intensity)))
+            end
         end
     end
     
@@ -442,11 +454,15 @@ function analyze_signal_quality(msi_data::MSIData; sample_indices::AbstractVecto
         return results
     end
     
+    results_lock = ReentrantLock()
+
     _iterate_spectra_fast(msi_data, sample_indices) do idx, mz, intensity
         if !isempty(intensity)
             # Noise estimation using MAD
             noise = mad(intensity, normalize=true)
-            push!(noise_levels, noise)
+            lock(results_lock) do
+                push!(noise_levels, noise)
+            end
             
             # FIX: More robust SNR calculation
             valid_intensity = intensity[intensity .> 0]  # Remove zeros
@@ -458,12 +474,16 @@ function analyze_signal_quality(msi_data::MSIData; sample_indices::AbstractVecto
                 if noise_robust > 0 && isfinite(signal_estimate)
                     snr_val = signal_estimate / noise_robust
                     # Cap unrealistic SNR values
-                    push!(snr_distribution, min(snr_val, 1e6))
+                    lock(results_lock) do
+                        push!(snr_distribution, min(snr_val, 1e6))
+                    end
                 end
             end
             
             # Total ion count
-            push!(tic_values, sum(intensity))
+            lock(results_lock) do
+                push!(tic_values, sum(intensity))
+            end
         end
     end
     
@@ -651,20 +671,23 @@ function analyze_peak_characteristics(msi_data::MSIData, instrument_analysis::Di
         peak_counts = Int[]
         fwhm_values = Float64[]
         
-        spectra_analyzed = 0
+        spectra_analyzed = Base.Threads.Atomic{Int}(0)
-        peaks_analyzed = 0
+        peaks_analyzed = Base.Threads.Atomic{Int}(0)
+        results_lock = ReentrantLock()
         
         _iterate_spectra_fast(msi_data, spectrum_indices) do idx, mz, intensity
             if length(mz) < 10  # Skip spectra with too few points
                 return
             end
             
-            spectra_analyzed += 1
+            Base.Threads.atomic_add!(spectra_analyzed, 1)
             meta = msi_data.spectra_metadata[idx]
             
             # Detect peaks with lower SNR threshold to find more peaks
             peaks = detect_peaks_profile_core(mz, intensity; snr_threshold=2.0)
-            push!(peak_counts, length(peaks))
+            lock(results_lock) do
+                push!(peak_counts, length(peaks))
+            end
             
             if !isempty(peaks)
                 # Analyze the strongest 3 peaks per spectrum
@@ -679,15 +702,17 @@ function analyze_peak_characteristics(msi_data::MSIData, instrument_analysis::Di
                         if !isnan(fwhm_delta_m) && fwhm_delta_m > 0.001 && fwhm_delta_m < 0.5  # Reasonable range in Da
                             fwhm_ppm = 1e6 * fwhm_delta_m / peak.mz
                             if 5.0 < fwhm_ppm < 500.0  # Reasonable ppm range
-                                push!(peak_widths_ppm, fwhm_ppm)
+                                lock(results_lock) do
-                                push!(fwhm_values, fwhm_delta_m)
+                                    push!(peak_widths_ppm, fwhm_ppm)
+                                    push!(fwhm_values, fwhm_delta_m)
+                                    
+                                    r2 = _fit_gaussian_and_r2(mz, intensity, peak_idx, 5)
+                                    push!(r_squared_values, r2)
+                                end
+                                Base.Threads.atomic_add!(peaks_analyzed, 1)
                                 
-                                r2 = _fit_gaussian_and_r2(mz, intensity, peak_idx, 5)
+                                if peaks_analyzed[] <= 3
-                                push!(r_squared_values, r2)
+                                    println("DEBUG: Peak at m/z $(peak.mz), FWHM = $(fwhm_ppm) ppm, R² = $(r_squared_values[end])")
-                                peaks_analyzed += 1
-                                
-                                if peaks_analyzed <= 3
-                                    println("DEBUG: Peak at m/z $(peak.mz), FWHM = $(fwhm_ppm) ppm, R² = $r2")
                                 end
                             end
                         end
@@ -698,7 +723,7 @@ function analyze_peak_characteristics(msi_data::MSIData, instrument_analysis::Di
             end
         end
         
-        println("DEBUG: Analyzed $peaks_analyzed peaks from $spectra_analyzed spectra")
+        println("DEBUG: Analyzed $(peaks_analyzed[]) peaks from $(spectra_analyzed[]) spectra")
         
         if !isempty(peak_widths_ppm)
             results[:mean_fwhm_ppm] = mean(peak_widths_ppm)
@@ -727,9 +752,13 @@ function analyze_peak_characteristics(msi_data::MSIData, instrument_analysis::Di
         
         peak_counts = Int[]
         
+        results_lock = ReentrantLock()
+
         _iterate_spectra_fast(msi_data, spectrum_indices) do idx, mz, intensity
             if !isempty(mz)
-                push!(peak_counts, length(mz))
+                lock(results_lock) do
+                    push!(peak_counts, length(mz))
+                end
             end
         end
         
@@ -1285,18 +1314,16 @@ function _fit_gaussian_and_r2(mz::AbstractVector{<:Real}, intensity::AbstractVec
     end_idx = min(n, peak_idx + half_window)
     
     # Ensure there's enough data to fit
-    if (end_idx - start_idx + 1) < 3
+    count = end_idx - start_idx + 1
+    if count < 3
         return 0.0
     end
 
-    x_data = mz[start_idx:end_idx]
-    y_data = intensity[start_idx:end_idx]
-
     # Estimate Gaussian parameters
     # Amplitude (A): peak intensity
-    A_est = intensity[peak_idx]
+    A_est = float(intensity[peak_idx])
     # Mean (μ): m/z at peak intensity
-    mu_est = mz[peak_idx]
+    mu_est = float(mz[peak_idx])
     # Standard deviation (σ): related to FWHM. FWHM = 2 * sqrt(2 * ln(2)) * σ ≈ 2.355 * σ
     # So, σ ≈ FWHM / 2.355
     fwhm_delta_m = _calculate_fwhm_delta_m(mz, intensity, peak_idx)
@@ -1310,19 +1337,27 @@ function _fit_gaussian_and_r2(mz::AbstractVector{<:Real}, intensity::AbstractVec
         return 0.0
     end
 
-    # Gaussian function
+    # Calculate SS_res, mean_y in a single pass to avoid allocations
-    gaussian(x, A, mu, sigma) = A * exp.(-(x .- mu).^2 ./ (2 * sigma^2))
+    SS_res = 0.0
+    sum_y = 0.0
+    
+    @inbounds for i in start_idx:end_idx
+        x_val = float(mz[i])
+        y_val = float(intensity[i])
+        
+        # Gaussian function estimate
+        y_est = A_est * exp(-((x_val - mu_est)^2) / (2 * sigma_est^2))
+        
+        SS_res += (y_val - y_est)^2
+        sum_y += y_val
+    end
 
-    # Generate estimated Gaussian curve
+    mean_y = sum_y / count
-    y_est = gaussian(x_data, A_est, mu_est, sigma_est)
+    SS_tot = 0.0
-
+    
-    # Calculate pseudo R-squared
+    @inbounds for i in start_idx:end_idx
-    # R^2 = 1 - (SS_res / SS_tot)
+        SS_tot += (float(intensity[i]) - mean_y)^2
-    # SS_res = sum((y_data - y_est).^2)
+    end
-    # SS_tot = sum((y_data - mean(y_data)).^2)
-
-    SS_res = sum((y_data .- y_est).^2)
-    SS_tot = sum((y_data .- mean(y_data)).^2)
 
     if SS_tot == 0
         return 1.0 # Perfect fit if all y_data are the same

src/Preprocessing.jl

@@ -12,6 +12,7 @@ generation.
 # =============================================================================
 
 using Statistics  # For mean, median
+using SparseArrays
 using StatsBase   # For mad (Median Absolute Deviation)
 using SavitzkyGolay # For SavitzkyGolay filtering
 using Dates         # For now()
@@ -40,7 +41,7 @@ A struct to hold the final feature matrix generated from the preprocessing pipel
 - `sample_ids::Vector{Int}`: A vector of identifiers for each sample (row) in the `matrix`.
 """
 struct FeatureMatrix
-    matrix::Array{Float64,2}
+    matrix::AbstractMatrix{Float64}
     mz_bins::Vector{Tuple{Float64,Float64}}
     sample_ids::Vector{Int}
 end
@@ -487,32 +488,23 @@ Estimates the baseline of a spectrum using the SNIP algorithm (internal implemen
 function _snip_baseline_impl(y::AbstractVector{<:Real}; iterations::Int=100)
     n = length(y)
     
-    # Initialize two buffers. b1 holds the current baseline estimate, b2 for the next.
+    # Initialize the baseline estimate array once
-    # Always convert to Float64 to ensure type stability and avoid copying if already correct type
     b1 = collect(float.(y))
-    b2 = similar(b1)        
     
-    current_b = b1
-    next_b = b2
-
     for k in 1:iterations
-        # Calculate next baseline estimate into `next_b` based on `current_b`
+        prev_val = b1[1]
-        # Boundary conditions
+        b1[1] = min(b1[1], b1[2])
-        if n > 1
-            next_b[1] = min(current_b[1], current_b[2])
-            next_b[n] = min(current_b[n], current_b[n-1])
-        end
-
-        @inbounds for i in 2:n-1
-            next_b[i] = min(current_b[i], 0.5 * (current_b[i-1] + current_b[i+1]))
-        end
         
-        # Swap references for the next iteration (no data copy here)
+        @inbounds for i in 2:n-1
-        current_b, next_b = next_b, current_b 
+            curr_val = b1[i]
+            b1[i] = min(curr_val, 0.5 * (prev_val + b1[i+1]))
+            prev_val = curr_val
+        end
+        b1[n] = min(b1[n], prev_val)
     end
     
-    # Return the final baseline estimate (which is in current_b after the last swap)
+    # Return the final baseline estimate
-    return current_b
+    return b1
 end
 
 """
@@ -720,18 +712,32 @@ function detect_peaks_profile_core(mz::AbstractVector{<:Real}, y::AbstractVector
     n = length(y)
     n < 3 && return NamedTuple{(:mz, :intensity, :fwhm, :shape_r2, :snr, :prominence), Tuple{Float64, Float64, Float64, Float64, Float64, Float64}}[]
     
-    noise_level = mad(y, normalize=true) + eps(Float64)
+    # Fast, non-allocating noise estimation
-    ys = smooth_spectrum_core(y; method=:savitzky_golay, window=max(5, 2*half_window+1), order=2) # Use smoothed data for detection
+    mean_y = sum(y) / n
+    noise_level = (sum(abs.(y .- mean_y)) / n) * 1.5 + eps(Float64)
+    
+    ys = smooth_spectrum_core(y; method=:savitzky_golay, window=max(5, 2*half_window+1), order=2)
 
     candidate_peak_indices = Int[]
-    for i in 2:n-1
+    sizehint!(candidate_peak_indices, div(n, 10)) # Pre-allocate memory capacity
+    
+    @inbounds for i in 2:n-1
         left = max(1, i - half_window)
         right = min(n, i + half_window)
         
-        # Prominence check
+        # Avoid @view allocation in tight loop by manually computing minimums and maximums
-        prominence = ys[i] - max(minimum(@view ys[left:i]), minimum(@view ys[i:right]))
+        min_left = ys[left]
+        for j in left:i; min_left = min(min_left, ys[j]); end
         
-        if ys[i] >= maximum(@view ys[left:right]) && 
+        min_right = ys[i]
+        for j in i:right; min_right = min(min_right, ys[j]); end
+        
+        prominence = ys[i] - max(min_left, min_right)
+        
+        max_local = ys[left]
+        for j in left:right; max_local = max(max_local, ys[j]); end
+        
+        if ys[i] >= max_local && 
            (ys[i] > snr_threshold * noise_level) &&
            (prominence > min_peak_prominence * ys[i])
             push!(candidate_peak_indices, i)
@@ -768,7 +774,14 @@ function detect_peaks_profile_core(mz::AbstractVector{<:Real}, y::AbstractVector
         
         left = max(1, p_idx - half_window)
         right = min(n, p_idx + half_window)
-        prominence = ys[p_idx] - max(minimum(@view ys[left:p_idx]), minimum(@view ys[p_idx:right]))
+        
+        min_left = ys[left]
+        for j in left:p_idx; min_left = min(min_left, ys[j]); end
+        
+        min_right = ys[p_idx]
+        for j in p_idx:right; min_right = min(min_right, ys[j]); end
+        
+        prominence = ys[p_idx] - max(min_left, min_right)
 
         push!(detected_peaks, (mz=peak_mz, intensity=peak_int, fwhm=fwhm_ppm, shape_r2=shape_r2, snr=peak_snr, prominence=prominence))
     end

src/PreprocessingPipeline.jl

@@ -33,7 +33,15 @@ function apply_baseline_correction(spectra::Vector{MutableSpectrum}, params::Dic
 
     Threads.@threads for s in spectra
         if validate_spectrum(s.mz, s.intensity)
+            original_length = length(s.intensity)
             baseline = apply_baseline_correction_core(s.intensity; method=method, iterations=iterations, window=window)
+            
+            # CRITICAL: Ensure baseline has the same length as intensity
+            if length(baseline) != original_length
+                @warn "Baseline correction: length mismatch for spectrum $(s.id). baseline=$(length(baseline)), intensity=$original_length. Skipping this spectrum."
+                continue
+            end
+            
             s.intensity = max.(0.0, s.intensity .- baseline)
         end
     end
@@ -53,7 +61,16 @@ function apply_intensity_transformation(spectra::Vector{MutableSpectrum}, params
 
     Threads.@threads for s in spectra
         if validate_spectrum(s.mz, s.intensity)
-            s.intensity = transform_intensity_core(s.intensity; method=method)
+            original_length = length(s.intensity)
+            transformed = transform_intensity_core(s.intensity; method=method)
+            
+            # CRITICAL: Ensure transformation preserves array length
+            if length(transformed) != original_length
+                @warn "Intensity transformation: length mismatch for spectrum $(s.id). transformed=$(length(transformed)), original=$original_length. Skipping this spectrum."
+                continue
+            end
+            
+            s.intensity = transformed
         end
     end
 end
@@ -76,8 +93,16 @@ function apply_smoothing(spectra::Vector{MutableSpectrum}, params::Dict)
 
     Threads.@threads for s in spectra
         if validate_spectrum(s.mz, s.intensity)
-            smoothed_intensity = max.(0.0, smooth_spectrum_core(s.intensity; method=method, window=window, order=order))
+            original_length = length(s.intensity)
-            s.intensity = smoothed_intensity
+            smoothed_intensity = smooth_spectrum_core(s.intensity; method=method, window=window, order=order)
+            
+            # CRITICAL: Ensure smoothing preserves array length  
+            if length(smoothed_intensity) != original_length
+                @warn "Smoothing: length mismatch for spectrum $(s.id). smoothed=$(length(smoothed_intensity)), original=$original_length, method=$method, window=$window. Skipping this spectrum."
+                continue
+            end
+            
+            s.intensity = max.(0.0, smoothed_intensity)
         end
     end
 end
@@ -104,6 +129,7 @@ function apply_peak_picking(spectra::Vector{MutableSpectrum}, params::Dict)
 
     Threads.@threads for s in spectra
         if validate_spectrum(s.mz, s.intensity)
+            # old_len = length(s.peaks)
             if method == :profile
                 s.peaks = detect_peaks_profile_core(s.mz, s.intensity; snr_threshold=snr_threshold, half_window=half_window, min_peak_prominence=min_peak_prominence, merge_peaks_tolerance=merge_peaks_tolerance)
             elseif method == :wavelet
@@ -113,6 +139,7 @@ function apply_peak_picking(spectra::Vector{MutableSpectrum}, params::Dict)
             else
                 s.peaks = detect_peaks_profile_core(s.mz, s.intensity; snr_threshold=snr_threshold, half_window=half_window)
             end
+            # @info "Spectrum $(s.id): detected $(length(s.peaks)) peaks"
         else
             s.peaks = []
         end
@@ -182,12 +209,21 @@ function apply_calibration(spectra::Vector{MutableSpectrum}, params::Dict, refer
     Threads.@threads for i in 1:length(spectra)
         s = spectra[i]
         if validate_spectrum(s.mz, s.intensity)
+            original_length = length(s.mz)
             matched_peaks = find_calibration_peaks_core(s.mz, s.intensity, reference_masses; ppm_tolerance=ppm_tolerance)
             if length(matched_peaks) >= 2
                 measured = sort(collect(values(matched_peaks)))
                 theoretical = sort(collect(keys(matched_peaks)))
                 itp = linear_interpolation(measured, theoretical, extrapolation_bc=Line())
-                s.mz = itp(s.mz) # Modify mz-axis in-place
+                new_mz = itp(s.mz)
+                
+                # CRITICAL: Ensure m/z axis preserves array length
+                if length(new_mz) != original_length
+                    @warn "Calibration: length mismatch for spectrum $(s.id). new_mz=$(length(new_mz)), original=$original_length. Skipping this spectrum."
+                    continue
+                end
+                
+                s.mz = new_mz # Modify mz-axis in-place
             else
                 @warn "Spectrum $(s.id): insufficient reference peaks ($(length(matched_peaks)) found), skipping calibration."
             end
@@ -233,7 +269,16 @@ function apply_peak_alignment(spectra::Vector{MutableSpectrum}, params::Dict)
         current_peaks_mz = [p.mz for p in s.peaks]
         alignment_func = align_peaks_lowess_core(ref_peaks_mz, current_peaks_mz; method=method, tolerance=tolerance, tolerance_unit=tolerance_unit)
         
-        s.mz = alignment_func.(s.mz) # Update m/z axis
+        original_length = length(s.mz)
+        new_mz = alignment_func.(s.mz)
+        
+        # CRITICAL: Ensure alignment preserves array length
+        if length(new_mz) != original_length
+            @warn "Peak alignment: m/z length mismatch for spectrum $(s.id). new_mz=$(length(new_mz)), original=$original_length. Skipping this spectrum."
+            continue
+        end
+        
+        s.mz = new_mz # Update m/z axis
         
         # Update peak m/z values
         for i in 1:length(s.peaks)
@@ -258,7 +303,16 @@ function apply_normalization(spectra::Vector{MutableSpectrum}, params::Dict)
 
     Threads.@threads for s in spectra
         if validate_spectrum(s.mz, s.intensity)
-            s.intensity = apply_normalization_core(s.intensity; method=method)
+            original_length = length(s.intensity)
+            normalized = apply_normalization_core(s.intensity; method=method)
+            
+            # CRITICAL: Ensure normalization preserves array length
+            if length(normalized) != original_length
+                @warn "Normalization: length mismatch for spectrum $(s.id). normalized=$(length(normalized)), original=$original_length. Skipping this spectrum."
+                continue
+            end
+            
+            s.intensity = normalized
         end
     end
 end

src/ResourcePool.jl

@@ -0,0 +1,79 @@
+# src/ResourcePool.jl
+using Base.Threads
+
+"""
+    ResourcePool{T}
+
+A thread-safe pool for reusing objects of type `T` to minimize allocations and GC pressure.
+Specifically designed for high-performance computing tasks where large buffers are needed
+repeatedly across multiple threads.
+
+# Fields:
+- `pool::Vector{T}`: The underlying storage for idle resources.
+- `lock::ReentrantLock`: Ensures thread-safe access to the pool.
+- `max_size::Int`: Maximum number of resources to hold in the pool.
+- `constructor::Function`: A function to create a new resource if the pool is empty.
+"""
+mutable struct ResourcePool{T}
+    pool::Vector{T}
+    lock::ReentrantLock
+    max_size::Int
+    constructor::Function
+end
+
+
+
+"""
+    ResourcePool{T}(constructor::Function; max_size::Int=2 * nthreads())
+
+Creates a new `ResourcePool` for resources of type `T`.
+"""
+function ResourcePool{T}(constructor::Function; max_size::Int=2 * nthreads()) where T
+    return ResourcePool{T}(T[], ReentrantLock(), max_size, constructor)
+end
+
+"""
+    acquire(pool::ResourcePool{T}) -> T
+
+Retrieves a resource from the pool. If the pool is empty, a new resource is created
+using the constructor.
+"""
+function acquire(pool::ResourcePool{T}) where T
+    lock(pool.lock) do
+        if !isempty(pool.pool)
+            return pop!(pool.pool)
+        end
+    end
+    # Create new resource outside of lock to minimize contention
+    return pool.constructor()
+end
+
+"""
+    release!(pool::ResourcePool{T}, resource::T)
+
+Returns a resource to the pool for later reuse. If the pool is already at `max_size`,
+the resource is allowed to be garbage collected.
+"""
+function release!(pool::ResourcePool{T}, resource::T) where T
+    lock(pool.lock) do
+        if length(pool.pool) < pool.max_size
+            push!(pool.pool, resource)
+        end
+    end
+    return nothing
+end
+
+"""
+    with_resource(f::Function, pool::ResourcePool{T})
+
+Acquires a resource from the pool, executes the function `f(resource)`, and
+automatically releases the resource back to the pool when finished.
+"""
+function with_resource(f::Function, pool::ResourcePool{T}) where T
+    resource = acquire(pool)
+    try
+        return f(resource)
+    finally
+        release!(pool, resource)
+    end
+end

src/StreamingKernels.jl

@@ -0,0 +1,317 @@
+# src/StreamingKernels.jl
+# ============================================================================
+# In-Place Spectral Kernels for the Streaming Pipeline
+# 
+# These functions operate on raw (mz, intensity) views from the Sprint 1
+# Mmap engine. They write results back to the input buffers using .= to
+# achieve zero-allocation processing per spectrum.
+#
+# Design contract:
+#   - All !-suffixed functions modify their arguments in-place
+#   - If a kernel needs temporary storage, it borrows from data.resource_pool
+#   - No function creates MutableSpectrum objects
+# ============================================================================
+
+using Statistics: mean, median
+
+# =============================================================================
+# Category A: Purely Streamable Kernels
+# =============================================================================
+
+"""
+    normalize_inplace!(intensity::AbstractVector{<:Real}, method::Symbol)
+
+Normalizes intensity values in-place. Supports :tic, :median, :rms.
+Zero-allocation for the normalization itself.
+"""
+@inline function normalize_inplace!(intensity::AbstractVector{<:Real}, method::Symbol)
+    if method === :tic
+        s = sum(intensity)
+        if s > 0
+            intensity ./= s
+        end
+    elseif method === :median
+        m = median(intensity)
+        if m > 0
+            intensity ./= m
+        end
+    elseif method === :rms
+        s = sqrt(sum(abs2, intensity) / length(intensity))
+        if s > 0
+            intensity ./= s
+        end
+    end
+    return intensity
+end
+
+"""
+    transform_inplace!(intensity::AbstractVector{Float64}, method::Symbol)
+
+Applies intensity transformation in-place. Supports :sqrt, :log1p, :log, :log2, :log10.
+"""
+@inline function transform_inplace!(intensity::AbstractVector{Float64}, method::Symbol)
+    if method === :sqrt
+        @inbounds @simd for i in eachindex(intensity)
+            intensity[i] = sqrt(max(0.0, intensity[i]))
+        end
+    elseif method === :log1p
+        @inbounds @simd for i in eachindex(intensity)
+            intensity[i] = log1p(max(0.0, intensity[i]))
+        end
+    elseif method === :log
+        @inbounds @simd for i in eachindex(intensity)
+            intensity[i] = log(max(eps(Float64), intensity[i]))
+        end
+    elseif method === :log2
+        @inbounds @simd for i in eachindex(intensity)
+            intensity[i] = log2(max(eps(Float64), intensity[i]))
+        end
+    elseif method === :log10
+        @inbounds @simd for i in eachindex(intensity)
+            intensity[i] = log10(max(eps(Float64), intensity[i]))
+        end
+    end
+    return intensity
+end
+
+"""
+    smooth_inplace!(intensity::AbstractVector{Float64}, data::MSIData;
+                    method::Symbol=:savitzky_golay, window::Int=9, order::Int=2)
+
+Smooths intensity in-place using a temporary buffer from the resource pool.
+The SavitzkyGolay library allocates internally, but we copy the result back
+to the original buffer and return the pool buffer.
+"""
+function smooth_inplace!(intensity::AbstractVector{Float64}, scratch::AbstractVector{Float64}, data::MSIData;
+                         method::Symbol=:savitzky_golay, window::Int=9, order::Int=2)
+    n = length(intensity)
+    if n < 3
+        return intensity
+    end
+
+    if method === :savitzky_golay
+        win = isodd(window) ? window : window + 1
+        if n < win
+            return intensity
+        end
+        # SavitzkyGolay handles its own math but causes mild allocation.
+        res = SavitzkyGolay.savitzky_golay(collect(intensity), win, order)
+        @inbounds for i in eachindex(intensity)
+            intensity[i] = max(0.0, res.y[i])
+        end
+    elseif method === :moving_average
+        copyto!(scratch, intensity)
+        
+        half_w = div(window, 2)
+        @inbounds for i in 1:n
+            s_idx = max(1, i - half_w)
+            e_idx = min(n, i + half_w)
+            s = 0.0
+            @simd for j in s_idx:e_idx
+                s += scratch[j]
+            end
+            intensity[i] = max(0.0, s / (e_idx - s_idx + 1))
+        end
+    end
+    return intensity
+end
+
+"""
+    baseline_subtract_inplace!(intensity::AbstractVector{Float64}, data::MSIData;
+                               method::Symbol=:snip, iterations::Int=100, window::Int=20)
+
+Subtracts baseline from intensity in-place. Uses two pool buffers for the
+SNIP ping-pong iteration to avoid any heap allocation in the hot loop.
+"""
+function baseline_subtract_inplace!(intensity::AbstractVector{Float64}, scratch::AbstractVector{Float64}, data::MSIData;
+                                    method::Symbol=:snip, iterations::Int=100, window::Int=20)
+    n = length(intensity)
+    if n < 3
+        return intensity
+    end
+    
+    if method === :snip
+        copyto!(scratch, intensity)
+        
+        for k in 1:iterations
+            prev_val = scratch[1]
+            scratch[1] = min(scratch[1], scratch[2])
+            
+            @inbounds for i in 2:n-1
+                curr_val = scratch[i]
+                scratch[i] = min(curr_val, 0.5 * (prev_val + scratch[i+1]))
+                prev_val = curr_val
+            end
+            scratch[n] = min(scratch[n], prev_val)
+        end
+        
+        @inbounds @simd for i in 1:n
+            intensity[i] = max(0.0, intensity[i] - scratch[i])
+        end
+    elseif method === :convex_hull
+        baseline = convex_hull_baseline(intensity)
+        @inbounds @simd for i in eachindex(intensity)
+            intensity[i] = max(0.0, intensity[i] - baseline[i])
+        end
+    elseif method === :median
+        baseline = median_baseline(intensity; window=window)
+        @inbounds @simd for i in eachindex(intensity)
+            intensity[i] = max(0.0, intensity[i] - baseline[i])
+        end
+    end
+    return intensity
+end
+
+"""
+    detect_peaks_streaming(mz::AbstractVector, intensity::AbstractVector;
+                          method::Symbol=:profile, snr_threshold::Float64=3.0,
+                          half_window::Int=10, min_peak_prominence::Float64=0.1,
+                          merge_peaks_tolerance::Float64=0.002)
+
+Detects peaks and returns a vector of (mz, intensity) tuples.
+This delegates to existing _core functions but returns a lightweight format
+suitable for sparse accumulation (no NamedTuple overhead in the hot path).
+"""
+function detect_peaks_streaming(callback::Function, mz::AbstractVector{Float64}, intensity::AbstractVector{Float64}, scratch::AbstractVector{Float64};
+                               method::Symbol=:profile, snr_threshold::Float64=3.0,
+                               half_window::Int=10, min_peak_prominence::Float64=0.1,
+                               merge_peaks_tolerance::Float64=0.002)
+    n = length(intensity)
+    if n < 3
+        return
+    end
+
+    if method === :profile || method === :wavelet
+        # Zero-allocation noisy estimation (using mean of bottom half)
+        sum_i = 0.0
+        @simd for i in 1:n
+            sum_i += intensity[i]
+        end
+        mean_i = sum_i / n
+        
+        sum_noise = 0.0
+        count_noise = 0
+        @inbounds for i in 1:n
+            if intensity[i] < mean_i
+                sum_noise += intensity[i]
+                count_noise += 1
+            end
+        end
+        # Use * 1.5 as an approximation to MAD
+        noise_level = count_noise > 0 ? (sum_noise / count_noise) * 1.5 + eps(Float64) : mean_i + eps(Float64)
+        
+        # We will use the scratch buffer to store candidate indices to avoid allocating `Int[]`
+        # Because scratch is Float64, we can safely store integer indices up to 2^53 exactly.
+        num_candidates = 0
+        
+        @inbounds for i in 2:n-1
+            if intensity[i] > snr_threshold * noise_level
+                left = max(1, i - half_window)
+                right = min(n, i + half_window)
+                
+                is_max = true
+                for j in left:right
+                    if intensity[j] > intensity[i]
+                        is_max = false
+                        break
+                    end
+                end
+                
+                if is_max
+                    # Compute prominence
+                    min_left = intensity[i]
+                    for j in left:i
+                        if intensity[j] < min_left
+                            min_left = intensity[j]
+                        end
+                    end
+                    min_right = intensity[i]
+                    for j in i:right
+                        if intensity[j] < min_right
+                            min_right = intensity[j]
+                        end
+                    end
+                    
+                    prominence = intensity[i] - max(min_left, min_right)
+                    
+                    if prominence > min_peak_prominence * intensity[i]
+                        num_candidates += 1
+                        scratch[num_candidates] = i
+                    end
+                end
+            end
+        end
+        
+        # Merge close peaks
+        if num_candidates > 0
+            if merge_peaks_tolerance > 0
+                last_idx = trunc(Int, scratch[1])
+                # We emit the first peak lazily down below, so let's compact them in place
+                num_merged = 1
+                
+                for i in 2:num_candidates
+                    idx = trunc(Int, scratch[i])
+                    if (mz[idx] - mz[last_idx]) > merge_peaks_tolerance
+                        num_merged += 1
+                        scratch[num_merged] = idx
+                        last_idx = idx
+                    elseif intensity[idx] > intensity[last_idx]
+                        scratch[num_merged] = idx
+                        last_idx = idx
+                    end
+                end
+                num_candidates = num_merged
+            end
+            
+            # Emit merged peaks
+            for i in 1:num_candidates
+                idx = trunc(Int, scratch[i])
+                callback(mz[idx], intensity[idx])
+            end
+        end
+        
+    elseif method === :centroid
+        # Just use any value over snr_threshold * mean_noise
+        sum_i = sum(intensity)
+        mean_i = sum_i / n
+        noise_level = mean_i + eps(Float64)
+        
+        @inbounds for i in 1:n
+            if intensity[i] > snr_threshold * noise_level
+                callback(mz[i], intensity[i])
+            end
+        end
+    end
+end
+
+# =============================================================================
+# Category B: Conditionally Streamable Kernels (Fixed-Reference)
+# =============================================================================
+
+"""
+    calibrate_inplace!(mz::Vector{Float64}, intensity::AbstractVector,
+                       reference_masses::Vector{Float64}; ppm_tolerance::Float64=20.0)
+
+Calibrates the m/z axis in-place using a fixed dictionary of internal standard
+reference masses. This is streamable because the reference is constant.
+
+Returns `true` if calibration was applied, `false` if insufficient peaks were found.
+"""
+function calibrate_inplace!(mz::Vector{Float64}, intensity::AbstractVector,
+                            reference_masses::Vector{Float64}; ppm_tolerance::Float64=20.0)
+    matched_peaks = find_calibration_peaks_core(mz, intensity, reference_masses; 
+                                                 ppm_tolerance=ppm_tolerance)
+    if length(matched_peaks) < 2
+        return false # Insufficient reference peaks
+    end
+    
+    measured = sort(collect(values(matched_peaks)))
+    theoretical = sort(collect(keys(matched_peaks)))
+    itp = linear_interpolation(measured, theoretical, extrapolation_bc=Line())
+    
+    # Apply calibration in-place
+    @inbounds for i in eachindex(mz)
+        mz[i] = itp(mz[i])
+    end
+    return true
+end

src/StreamingPipeline.jl

@@ -0,0 +1,402 @@
+# src/StreamingPipeline.jl
+# ============================================================================
+# The Streaming Pipeline Executor
+#
+# This module provides `process_dataset!`, the Sprint 2 master function that
+# streams spectral data through an in-place kernel chain and accumulates
+# results into a SparseMatrixCSC without ever holding more than 1 spectrum
+# per thread in RAM.
+#
+# Architecture:
+#   1. _iterate_spectra_fast → Mmap zero-copy views
+#   2. copyto!(writable_buf, view) → makes mutable copy for kernels
+#   3. Kernel chain: smooth! → baseline! → peaks → bin
+#   4. Thread-local (I, J, V) sparse accumulators
+#   5. Final sparse(I, J, V, num_bins, num_spectra) assembly
+#
+# This works alongside the existing execute_full_preprocessing in
+# PreprocessingPipeline.jl — it does NOT replace the app.jl integration.
+# ============================================================================
+
+using SparseArrays
+using Printf
+
+# =============================================================================
+# Configuration Structs
+# =============================================================================
+
+"""
+    StreamingStep
+
+Represents a single step in the streaming pipeline.
+"""
+struct StreamingStep
+    name::Symbol
+    params::Dict{Symbol, Any}
+end
+
+"""
+    PipelineConfig
+
+Holds the complete configuration for a streaming pipeline execution.
+
+# Fields
+- `steps::Vector{StreamingStep}` — ordered sequence of processing steps
+- `reference_peaks::Vector{Float64}` — fixed m/z values for calibration (Category B)
+- `num_bins::Int` — number of bins for the output feature matrix
+- `min_peaks_per_bin::Int` — minimum peak count to keep a bin
+- `frequency_threshold::Float64` — minimum fraction of spectra a bin must appear in (0.0-1.0)
+
+# Example
+```julia
+config = PipelineConfig(
+    steps = [
+        StreamingStep(:smoothing, Dict(:method => :savitzky_golay, :window => 9, :order => 2)),
+        StreamingStep(:baseline_correction, Dict(:method => :snip, :iterations => 100)),
+        StreamingStep(:normalization, Dict(:method => :tic)),
+        StreamingStep(:peak_picking, Dict(:method => :profile, :snr_threshold => 3.0)),
+    ],
+    num_bins = 2000
+)
+```
+"""
+struct PipelineConfig
+    steps::Vector{StreamingStep}
+    reference_peaks::Vector{Float64}
+    num_bins::Int
+    min_peaks_per_bin::Int
+    frequency_threshold::Float64
+end
+
+# Convenience constructor with defaults
+function PipelineConfig(; steps::Vector{StreamingStep}=StreamingStep[],
+                         reference_peaks::Vector{Float64}=Float64[],
+                         num_bins::Int=2000,
+                         min_peaks_per_bin::Int=3,
+                         frequency_threshold::Float64=0.0)
+    return PipelineConfig(steps, reference_peaks, num_bins, min_peaks_per_bin, frequency_threshold)
+end
+
+# =============================================================================
+# Sparse Accumulator (Thread-Local)
+# =============================================================================
+
+"""
+    SparseAccumulator
+
+Thread-local accumulator for sparse matrix construction.
+Collects (row, col, val) triplets that will be assembled into
+a SparseMatrixCSC at the end of the pipeline.
+"""
+mutable struct SparseAccumulator
+    I::Vector{Int}    # Row indices (bin indices)
+    J::Vector{Int}    # Column indices (spectrum indices)
+    V::Vector{Float64} # Values (intensities)
+    lck::Base.Threads.SpinLock
+
+    function SparseAccumulator(capacity_hint::Int=10000)
+        acc = new(
+            Vector{Int}(undef, 0),
+            Vector{Int}(undef, 0),
+            Vector{Float64}(undef, 0),
+            Base.Threads.SpinLock()
+        )
+        sizehint!(acc.I, capacity_hint)
+        sizehint!(acc.J, capacity_hint)
+        sizehint!(acc.V, capacity_hint)
+        return acc
+    end
+end
+
+"""
+    accumulate!(acc::SparseAccumulator, spectrum_idx::Int, bin_indices::AbstractVector{Int},
+                intensities::AbstractVector{Float64})
+
+Appends peak data for one spectrum into the sparse accumulator.
+"""
+@inline function accumulate!(acc::SparseAccumulator, spectrum_idx::Int,
+                             bin_indices::AbstractVector{Int},
+                             intensities::AbstractVector{Float64})
+    n = length(bin_indices)
+    for k in 1:n
+        @inbounds begin
+            push!(acc.I, bin_indices[k])
+            push!(acc.J, spectrum_idx)
+            push!(acc.V, intensities[k])
+        end
+    end
+end
+
+# =============================================================================
+# The Pipeline Executor
+# =============================================================================
+
+"""
+    process_dataset!(data::MSIData, config::PipelineConfig;
+                     progress_callback::Union{Function, Nothing}=nothing,
+                     masked_indices::Union{AbstractVector{Int}, Nothing}=nothing)
+
+The Sprint 2 master streaming function. Processes an entire MSI dataset through
+a kernel chain without holding more than 1 spectrum per thread in RAM.
+
+# Returns
+- `SparseMatrixCSC{Float64, Int}`: The feature matrix (bins × spectra)
+- `Vector{Float64}`: The m/z bin centers
+
+# Architecture
+1. Ensures analytics are computed (for global m/z range)
+2. Creates thread-local SparseAccumulators
+3. Streams spectra via `_iterate_spectra_fast`
+4. Per spectrum: copy view → kernel chain → peak detect → bin → accumulate
+5. Merges accumulators → `sparse(I, J, V)`
+"""
+function process_dataset!(data::MSIData, config::PipelineConfig;
+                          progress_callback::Union{Function, Nothing}=nothing,
+                          masked_indices::Union{AbstractVector{Int}, Nothing}=nothing)
+    
+    # --- Step 1: Ensure analytics are computed (provides global m/z range) ---
+    if !is_set(data.analytics_ready)
+        println("Pre-computing analytics for streaming pipeline...")
+        precompute_analytics(data)
+    end
+    
+    # Determine global m/z range for binning
+    global_min_mz = Base.Threads.atomic_add!(data.global_min_mz, 0.0)
+    global_max_mz = Base.Threads.atomic_add!(data.global_max_mz, 0.0)
+    
+    if !isfinite(global_min_mz) || !isfinite(global_max_mz) || global_min_mz >= global_max_mz
+        @warn "Invalid global m/z range: [$global_min_mz, $global_max_mz]. Cannot bin peaks."
+        return spzeros(0, 0), Float64[]
+    end
+    
+    num_bins = config.num_bins
+    bin_edges = range(global_min_mz, stop=global_max_mz, length=num_bins + 1)
+    bin_centers = [(bin_edges[i] + bin_edges[i+1]) / 2 for i in 1:num_bins]
+    inv_bin_width = 1.0 / step(bin_edges)
+    
+    num_spectra = length(data.spectra_metadata)
+    indices_to_process = masked_indices === nothing ? nothing : masked_indices
+    
+    # --- Step 2: Create thread-local accumulators ---
+    n_threads = Base.Threads.nthreads()
+    accumulators = [SparseAccumulator(num_spectra * 10) for _ in 1:n_threads]
+    spectra_processed = Base.Threads.Atomic{Int}(0)
+    
+    # NEW: Create dedicated workspace buffers for each thread.
+    # This completely eliminates the need for acquire/release and prevents deadlocks.
+    workspaces_mz = [Vector{Float64}(undef, 0) for _ in 1:n_threads]
+    workspaces_int = [Vector{Float64}(undef, 0) for _ in 1:n_threads]
+    workspaces_scratch = [Vector{Float64}(undef, 0) for _ in 1:n_threads]
+    
+    # Pre-parse step configuration for fast dispatch in the hot loop
+    has_smoothing = false
+    has_baseline = false
+    has_normalization = false
+    has_transform = false
+    has_peak_picking = false
+    has_calibration = false
+    
+    smooth_params = Dict{Symbol, Any}()
+    baseline_params = Dict{Symbol, Any}()
+    norm_params = Dict{Symbol, Any}()
+    transform_params = Dict{Symbol, Any}()
+    peak_params = Dict{Symbol, Any}()
+    
+    for s in config.steps
+        if s.name === :smoothing
+            has_smoothing = true
+            smooth_params = s.params
+        elseif s.name === :baseline_correction
+            has_baseline = true
+            baseline_params = s.params
+        elseif s.name === :normalization
+            has_normalization = true
+            norm_params = s.params
+        elseif s.name === :stabilization || s.name === :intensity_transformation
+            has_transform = true
+            transform_params = s.params
+        elseif s.name === :peak_picking
+            has_peak_picking = true
+            peak_params = s.params
+        elseif s.name === :calibration
+            has_calibration = true
+        end
+    end
+    
+    reference_masses = config.reference_peaks
+    
+    # --- Step 3: Stream and process ---
+    start_time = time_ns()
+    
+    # Use let block to capture all variables cleanly for the closure
+    let data=data, accumulators=accumulators, spectra_processed=spectra_processed,
+        bin_edges=bin_edges, num_bins=num_bins, inv_bin_width=inv_bin_width,
+        global_min_mz=global_min_mz,
+        workspaces_mz=workspaces_mz, workspaces_int=workspaces_int, workspaces_scratch=workspaces_scratch,
+        has_smoothing=has_smoothing, has_baseline=has_baseline,
+        has_normalization=has_normalization, has_transform=has_transform,
+        has_peak_picking=has_peak_picking, has_calibration=has_calibration,
+        smooth_params=smooth_params, baseline_params=baseline_params,
+        norm_params=norm_params, transform_params=transform_params,
+        peak_params=peak_params, reference_masses=reference_masses
+        
+        _iterate_spectra_fast(data, indices_to_process) do idx, mz_view, int_view
+            thread_id = Base.Threads.threadid()
+            acc = accumulators[thread_id]
+            
+            # --- Grab Thread-Local Workspaces ---
+            # No locking, no blocking, guaranteed to be available
+            mz_buf = workspaces_mz[thread_id]
+            int_buf = workspaces_int[thread_id]
+            scratch_buf = workspaces_scratch[thread_id]
+            
+            resize!(mz_buf, length(mz_view))
+            resize!(int_buf, length(int_view))
+            resize!(scratch_buf, length(int_view))
+            copyto!(mz_buf, mz_view)
+            copyto!(int_buf, int_view)
+            
+            # --- Kernel Chain (in pipeline order) ---
+                
+                # Category B: Fixed-reference calibration
+                if has_calibration && !isempty(reference_masses)
+                    calibrate_inplace!(mz_buf, int_buf, reference_masses)
+                end
+                
+                # Category A: Intensity transformation
+                if has_transform
+                    transform_inplace!(int_buf, get(transform_params, :method, :sqrt))
+                end
+                
+                # Category A: Smoothing
+                if has_smoothing
+                    smooth_inplace!(int_buf, scratch_buf, data;
+                        method=get(smooth_params, :method, :savitzky_golay),
+                        window=get(smooth_params, :window, 9),
+                        order=get(smooth_params, :order, 2))
+                end
+                
+                # Category A: Baseline correction
+                if has_baseline
+                    baseline_subtract_inplace!(int_buf, scratch_buf, data;
+                        method=get(baseline_params, :method, :snip),
+                        iterations=get(baseline_params, :iterations, 100),
+                        window=get(baseline_params, :window, 20))
+                end
+                
+                # Category A: Normalization
+                if has_normalization
+                    normalize_inplace!(int_buf, get(norm_params, :method, :tic))
+                end
+                
+                # --- Peak Detection & Binning ---
+                if has_peak_picking
+                    detect_peaks_streaming(mz_buf, int_buf, scratch_buf;
+                        method=get(peak_params, :method, :profile),
+                        snr_threshold=Float64(get(peak_params, :snr_threshold, 3.0)),
+                        half_window=Int(get(peak_params, :half_window, 10)),
+                        min_peak_prominence=Float64(get(peak_params, :min_peak_prominence, 0.1)),
+                        merge_peaks_tolerance=Float64(get(peak_params, :merge_peaks_tolerance, 0.002))) do peak_mz, peak_int
+                        
+                        # Bin each discovered peak directly
+                        bin_idx = trunc(Int, (peak_mz - global_min_mz) * inv_bin_width) + 1
+                        bin_idx = clamp(bin_idx, 1, num_bins)
+                        
+                        push!(acc.I, bin_idx)
+                        push!(acc.J, idx)
+                        push!(acc.V, peak_int)
+                    end
+                else
+                    # No peak picking: bin raw intensity directly
+                    @inbounds for i in eachindex(mz_buf)
+                        bin_idx = trunc(Int, (mz_buf[i] - global_min_mz) * inv_bin_width) + 1
+                        bin_idx = clamp(bin_idx, 1, num_bins)
+                        
+                        push!(acc.I, bin_idx)
+                        push!(acc.J, idx)
+                        push!(acc.V, int_buf[i])
+                    end
+                end
+                
+                Base.Threads.atomic_add!(spectra_processed, 1)
+        end
+    end
+    
+    # --- Step 4: Merge thread-local accumulators ---
+    total_entries = sum(length(acc.I) for acc in accumulators)
+    merged_I = Vector{Int}(undef, total_entries)
+    merged_J = Vector{Int}(undef, total_entries)
+    merged_V = Vector{Float64}(undef, total_entries)
+    
+    offset = 0
+    for acc in accumulators
+        n = length(acc.I)
+        if n > 0
+            copyto!(merged_I, offset + 1, acc.I, 1, n)
+            copyto!(merged_J, offset + 1, acc.J, 1, n)
+            copyto!(merged_V, offset + 1, acc.V, 1, n)
+            offset += n
+        end
+    end
+    
+    # --- Step 5: Assemble sparse matrix ---
+    # Use max combiner: when multiple peaks map to the same bin for same spectrum,
+    # keep the maximum intensity
+    feature_matrix = sparse(merged_I, merged_J, merged_V, num_bins, num_spectra, max)
+    
+    # --- Step 6: Apply frequency threshold if configured ---
+    if config.frequency_threshold > 0.0
+        # Count how many spectra have a non-zero value in each bin
+        bin_presence = vec(sum(feature_matrix .> 0, dims=2))
+        min_count = ceil(Int, config.frequency_threshold * num_spectra)
+        keep_bins = findall(bin_presence .>= min_count)
+        feature_matrix = feature_matrix[keep_bins, :]
+        bin_centers = bin_centers[keep_bins]
+    end
+    
+    duration = (time_ns() - start_time) / 1e9
+    n_processed = spectra_processed[]
+    n_nonzeros = nnz(feature_matrix)
+    sparsity = 1.0 - n_nonzeros / (size(feature_matrix, 1) * size(feature_matrix, 2) + 1)
+    
+    @printf "Streaming pipeline complete: %d spectra processed in %.2f seconds.\n" n_processed duration
+    @printf "Feature matrix: %d bins × %d spectra, %d non-zeros (%.1f%% sparse)\n" size(feature_matrix, 1) size(feature_matrix, 2) n_nonzeros sparsity * 100
+    @printf "RAM: %.1f MB (vs %.1f MB dense)\n" (n_nonzeros * 16) / 1e6 (size(feature_matrix, 1) * size(feature_matrix, 2) * 8) / 1e6
+    
+    if progress_callback !== nothing
+        progress_callback(1.0)
+    end
+    
+    return feature_matrix, collect(Float64, bin_centers)
+end
+
+"""
+    save_sparse_matrix(matrix::SparseMatrixCSC, output_path::String)
+
+Exports a highly optimized SparseMatrixCSC array to disk using the standard
+Matrix Market Coordinate format (`.mtx`), guaranteeing no bottleneck or OOM crashes 
+for extremely large MS dataset persistence.
+"""
+function save_sparse_matrix(matrix::SparseMatrixCSC{Float64, Int}, output_path::String)
+    m, n = size(matrix)
+    nnz_val = nnz(matrix)
+    # Use streaming I/O with a large buffer for ultra-fast persistence
+    open(output_path, "w") do io
+        # Write Matrix Market Header
+        write(io, "%%MatrixMarket matrix coordinate real general\n")
+        write(io, "$m $n $nnz_val\n")
+        
+        # Directly extract CSC properties (O(1) memory, zero allocation)
+        row_indices = rowvals(matrix)
+        values_array = nonzeros(matrix)
+        
+        @inbounds for filter_j in 1:n
+            # nzrange returns the index bounds for non-zero elements in column 'j'
+            for idx in nzrange(matrix, filter_j)
+                i = row_indices[idx]
+                v = values_array[idx]
+                write(io, "$i $filter_j $v\n")
+            end
+        end
+    end
+end

src/mzML.jl

@@ -19,6 +19,24 @@ const DATA_FORMAT_ACCESSIONS = Dict{String, DataType}(
     "MS:1000523" => Float64
 )
 
+"""
+    parse_instrument_metadata_mzml(stream::IO)
+
+Parses instrument metadata from an mzML file.
+
+# Arguments
+
+- `stream::IO`: The stream to parse the metadata from.
+
+# Returns
+- `InstrumentMetadata`: The instrument metadata.
+
+# Example
+
+```julia
+parse_instrument_metadata_mzml(open("test.mzML"))
+```
+"""
 function parse_instrument_metadata_mzml(stream::IO)
     println("DEBUG: Starting mzML instrument metadata parsing...")
     # Initialize with default values from the InstrumentMetadata constructor
@@ -213,7 +231,7 @@ function get_spectrum_asset_metadata(stream::IO)
     #println("DEBUG: Exiting get_spectrum_asset_metadata.")
 
     # Create SpectrumAsset directly from the variables
-    return SpectrumAsset(data_format, compression_flag, binary_offset, encoded_length, axis)
+    return SpectrumAsset(data_format, compression_flag, binary_offset, encoded_length, axis, 0.0, 0.0)
 end
 
 # This function is updated to return the generic SpectrumMetadata struct
@@ -328,6 +346,19 @@ end
 
 Finds the index offset in an mzML file by reading from the end.
 Optimized version with better buffer management.
+
+# Arguments
+
+- `stream::IO`: The stream to parse the metadata from.
+
+# Returns
+- `Int64`: The index offset.
+
+# Example
+
+```julia
+find_index_offset(open("test.mzML"))
+```
 """
 function find_index_offset(stream::IO)::Int64
     file_size = filesize(stream)
@@ -363,38 +394,48 @@ then parses the metadata for each spectrum without loading the binary data.
 """
 function load_mzml_lazy(file_path::String; cache_size::Int=100)
     println("DEBUG: Opening file stream for $file_path")
-    ts_stream = ThreadSafeFileHandle(file_path, "r")
+    
+    # --- Handle Pool Optimization ---
+    # Open multiple handles to the .mzML file to avoid lock contention
+    num_handles = Threads.nthreads()
+    mzml_handles = [open(file_path, "r") for _ in 1:num_handles]
+    
+    # Use the first handle for initial parsing
+    primary_handle = mzml_handles[1]
 
     try
         # --- NEW: Parse instrument metadata from header ---
         println("DEBUG: Parsing instrument metadata from header...")
-        instrument_meta = parse_instrument_metadata_mzml(ts_stream.handle)
+        instrument_meta = parse_instrument_metadata_mzml(primary_handle)
 
-        println("--- Extracted Instrument Metadata ---")
+        seekstart(primary_handle) # Reset stream after header parsing
-        println("Resolution: ", instrument_meta.resolution)
-        println("Acquisition Mode (pre-check): ", instrument_meta.acquisition_mode)
-        println("Calibration Status: ", instrument_meta.calibration_status)
-        println("Instrument Model: ", instrument_meta.instrument_model)
-        println("Mass Accuracy (ppm): ", instrument_meta.mass_accuracy_ppm)
-        println("Laser Settings: ", instrument_meta.laser_settings)
-        println("Polarity: ", instrument_meta.polarity)
-        println("------------------------------------")
-
-        seekstart(ts_stream.handle) # Reset stream after header parsing
 
         println("DEBUG: Finding index offset...")
-        index_offset = find_index_offset(ts_stream.handle)
+        index_offset = find_index_offset(primary_handle)
+        
+        # --- NEW: Mmap Optimization with RAM Safety ---
+        mmap_data = nothing
+        try
+            file_size = filesize(file_path)
+
+            println("DEBUG: Memory mapping .mzML file...")
+            seekstart(primary_handle) # Anchor Mmap to the beginning of the file to prevent overflow
+            mmap_data = Mmap.mmap(primary_handle, Vector{UInt8}, (file_size,))
+            println("DEBUG: .mzML file mmapped successfully.")
+        catch e
+            @warn "Memory mapping failed for mzML, falling back to standard I/O: $e"
+        end
+
         println("DEBUG: Seeking to index list at offset $index_offset.")
-        seek(ts_stream.handle, index_offset)
+        seek(primary_handle, index_offset)
         
         println("DEBUG: Searching for '<index name=\"spectrum\">'.")
-        if find_tag(ts_stream.handle, r"<index\s+name=\"spectrum\"") === nothing
+        if find_tag(primary_handle, r"<index\s+name=\"spectrum\"") === nothing
             throw(FileFormatError("Could not find spectrum index."))
         end
-        println("DEBUG: Found spectrum index tag.")
 
         println("DEBUG: Parsing spectrum offsets...")
-        spectrum_offsets = parse_offset_list(ts_stream.handle)
+        spectrum_offsets = parse_offset_list(primary_handle)
         if isempty(spectrum_offsets)
             throw(FileFormatError("No spectrum offsets found."))
         end
@@ -402,31 +443,25 @@ function load_mzml_lazy(file_path::String; cache_size::Int=100)
         println("DEBUG: Found $num_spectra spectrum offsets.")
 
         println("DEBUG: Parsing metadata for each spectrum...")
-        # Pre-allocate the metadata vector for better performance
         spectra_metadata = Vector{SpectrumMetadata}(undef, num_spectra)
         
-        # Use @inbounds for faster indexing in the loop
         @inbounds for i in 1:num_spectra
-            spectra_metadata[i] = parse_spectrum_metadata(ts_stream.handle, spectrum_offsets[i])
+            spectra_metadata[i] = parse_spectrum_metadata(primary_handle, spectrum_offsets[i])
             
-            # Progress reporting for large files
             if i % 1000 == 0
                 println("DEBUG: Processed $i/$num_spectra spectra")
             end
         end
         println("DEBUG: Metadata parsing complete for all $num_spectra spectra.")
 
-        # Assuming uniform data formats, take from the first spectrum
+        # Inferred global formats from first spectrum
         first_meta = spectra_metadata[1]
         mz_format = first_meta.mz_asset.format
         intensity_format = first_meta.int_asset.format
-        println("DEBUG: Inferred global m/z format: $mz_format")
-        println("DEBUG: Inferred global intensity format: $intensity_format")
         
-        # --- NEW: Determine overall acquisition mode ---
+        # Determine overall acquisition mode ...
-        modes = [meta.mode for meta in spectra_metadata]
+        num_centroid = count(m -> m.mode == CENTROID, spectra_metadata)
-        num_centroid = count(m -> m == CENTROID, modes)
+        num_profile = count(m -> m.mode == PROFILE, spectra_metadata)
-        num_profile = count(m -> m == PROFILE, modes)
         
         acq_mode_symbol = if num_centroid > 0 && num_profile == 0
             :centroid
@@ -437,26 +472,28 @@ function load_mzml_lazy(file_path::String; cache_size::Int=100)
         else
             :unknown
         end
-        println("DEBUG: Inferred overall acquisition mode: $acq_mode_symbol (Centroid: $num_centroid, Profile: $num_profile)")
 
         final_instrument_meta = InstrumentMetadata(
             instrument_meta.resolution,
-            acq_mode_symbol, # Update with parsed mode
+            acq_mode_symbol, 
             instrument_meta.mz_axis_type,
             instrument_meta.calibration_status,
             instrument_meta.instrument_model,
             instrument_meta.mass_accuracy_ppm,
             instrument_meta.laser_settings,
             instrument_meta.polarity,
-            instrument_meta.vendor_preprocessing_steps # Add this new field
+            instrument_meta.vendor_preprocessing_steps
         )
 
-        source = MzMLSource(ts_stream, mz_format, intensity_format)
+        source = MzMLSource(mzml_handles, mz_format, intensity_format, mmap_data)
         println("DEBUG: Creating MSIData object.")
         return MSIData(source, spectra_metadata, final_instrument_meta, (0, 0), nothing, cache_size)
 
     catch e
-        close(ts_stream) # Ensure stream is closed on error
+        # Close all handles in the pool if initialization fails
+        for h in mzml_handles
+            isopen(h) && close(h)
+        end
         rethrow(e)
     end
 end

start_MSI_GUI.jl

@@ -6,21 +6,65 @@ ENV["GENIE_ENV"] = "dev"
 
 manifest_path = joinpath(@__DIR__, "Manifest.toml")
 
-# Only instantiate in development mode
+# Selective instantiation for faster startup
-if get(ENV, "GENIE_ENV", "dev") != "prod" || !isfile(manifest_path)
+if get(ENV, "GENIE_ENV", "dev") != "prod" && !isfile(manifest_path)
-    @info "Development environment detected. Instantiating packages..."
+    @info "Development environment detected and Manifest.toml missing. Instantiating packages..."
-    
-    if !isfile(manifest_path)
-        @info "Manifest.toml not found. Generating it based on Project.toml..."
-    end
-
     Pkg.resolve()
     Pkg.instantiate()
     Pkg.gc()
+elseif get(ENV, "GENIE_ENV", "dev") != "prod"
+    @info "Manifest.toml found. Skipping Pkg.instantiate() for faster boot. Delete Manifest.toml if you need to re-instantiate."
+    ENV["GENIE_ENV"] = "prod"
 end
 
 using Genie
 
+# --- Cross-Platform Startup Cleanup ---
+# Remove orphaned GenieSessionFileSession directories from previous runs.
+# These accumulate in the OS temp directory as jl_XXXXXX folders containing
+# serialized session files (64-char hex filenames). Over long sessions or
+# after crashes, they can consume gigabytes of disk space.
+function cleanup_orphaned_sessions()
+    tmp = Base.tempdir()
+    cleaned_count = 0
+    cleaned_bytes = 0
+
+    for entry in readdir(tmp; join=false)
+        # Only target directories matching Julia's temp naming pattern
+        startswith(entry, "jl_") || continue
+        full_path = joinpath(tmp, entry)
+        isdir(full_path) || continue
+
+        # Validate: a Genie session dir contains files with 64-char hex names
+        try
+            contents = readdir(full_path)
+            isempty(contents) && continue
+
+            # Check if at least one file matches the 64-char hex session ID pattern
+            is_session_dir = any(contents) do f
+                length(f) == 64 && all(c -> c in "0123456789abcdef", f)
+            end
+            is_session_dir || continue
+
+            # Safe to remove — this is an orphaned Genie session directory
+            dir_size = sum(filesize(joinpath(full_path, f)) for f in contents; init=0)
+            rm(full_path; recursive=true, force=true)
+            cleaned_count += 1
+            cleaned_bytes += dir_size
+        catch e
+            @debug "Skipping $entry during cleanup: $e"
+        end
+    end
+
+    if cleaned_count > 0
+        size_mb = round(cleaned_bytes / (1024^2), digits=1)
+        @info "Startup cleanup: removed $cleaned_count orphaned session dir(s), freed $(size_mb) MB"
+    end
+end
+
+cleanup_orphaned_sessions()
+
+
 # Load and configure Genie
 Genie.loadapp()
 

test/benchmark.jl

@@ -49,12 +49,16 @@ const BENCHMARK_CASES = [
     # BenchmarkCase("/path/to/your/small_file.imzML", 309.06, 0.1),
     # BenchmarkCase("/path/to/your/medium_file.imzML", 896.0, 1.0),
     # BenchmarkCase("/path/to/your/large_file.imzML", 100.0, 0.1),
-    BenchmarkCase("/home/pixel/Documents/Cinvestav_2025/Analisis/imzML_LA-ESI/180817_NEG_Thaliana_Leaf_bottom_1_0841.imzML",116.07,0.1, "Thaliana Leaf"),
+    # BenchmarkCase("/home/pixel/Documents/Cinvestav_2025/Analisis/imzML_LA-ESI/180817_NEG_Thaliana_Leaf_bottom_1_0841.imzML",116.07,0.1, "Thaliana Leaf"),
-    BenchmarkCase("/home/pixel/Documents/Cinvestav_2025/Analisis/imzML_LTP/ltpmsi-chilli.imzML",420,0.1, "Chilli Pepper"), # Chilli
+    BenchmarkCase("/home/pixel/Documents/Cinvestav_2025/Analisis/Chilli/ltpmsi-chilli.imzML",420,0.1, "Chilli Pepper"), # Chilli
-    BenchmarkCase("/home/pixel/Documents/Cinvestav_2025/Analisis/imzML_DESI/ColAd_Individual/40TopL,10TopR,30BottomL,20BottomR/40TopL,10TopR,30BottomL,20BottomR-centroid.imzML",885.5,0.1, "Colon Cancer Human"), # Human Cancer
+    BenchmarkCase("/home/pixel/Documents/Cinvestav_2025/Analisis/imzML_DESI/40TopL,10TopR,30BottomL,20BottomR/40TopL,10TopR,30BottomL,20BottomR-centroid.imzML",885.5,0.1, "Colon Cancer Human"), # Human Cancer
     BenchmarkCase("/home/pixel/Documents/Cinvestav_2025/Analisis/imzML_AP_SMALDI/HR2MSImouseurinarybladderS096.imzML", 716.053,0.1, "Mouse Urinary Bladder"), # Mouse bladder
+    BenchmarkCase("/home/pixel/Documents/Cinvestav_2025/Analisis/Leafs/CE1_Leaf_R3.imzML",306.1,0.1, "Leaf"),
     BenchmarkCase("/home/pixel/Documents/Cinvestav_2025/Analisis/Liv2_imzML_TIMSConvert-selected/Liv2.imzML",796.18,0.1, "Liver Cut"), #Lib2
-    BenchmarkCase("/home/pixel/Documents/Cinvestav_2025/Analisis/salida/Stomach_DHB_uncompressed.imzML",804.3,0.1, "Mouse Stomach"), # Mouse Stomach
+    BenchmarkCase("/home/pixel/Documents/Cinvestav_2025/Analisis/salida/Stomach_DHB_uncompressed.imzML",804.3,0.1, "Mouse Stomach 2GB"), # Mouse Stomach
+    BenchmarkCase("/home/pixel/Documents/Cinvestav_2025/Analisis/salida_Tims/Stomach_DHB.imzML",804.3,0.1, "Mouse Stomach 4GB"), # Mouse Stomach
+    
+
 ]
 
 const NUM_REPETITIONS = 50 # Number of times to generate the image for averaging

test/benchmark_v3.jl

@@ -0,0 +1,215 @@
+# test/benchmark_v3.jl
+
+# ===================================================================
+# High-Precision Performance Benchmark Suite (v3)
+# ===================================================================
+# This script integrates the robust statistical sampling of `BenchmarkTools`
+# with the visual comparative mechanics against the legacy library.
+# It captures the paradigm shift from "Time per slice" to "Pipeline Throughput".
+# ===================================================================
+
+using Pkg
+Pkg.activate(joinpath(@__DIR__, ".."))
+
+using BenchmarkTools
+using DataFrames
+using CSV
+using CairoMakie
+using Statistics
+using MSI_src
+using julia_mzML_imzML
+
+struct BenchmarkCase
+    filepath::String
+    mz_value::Float64
+    mz_tolerance::Float64
+    name::String
+end
+
+const RESULTS_DIR = joinpath(@__DIR__, "results")
+
+# List to benchmark
+const BENCHMARK_CASES = [
+    BenchmarkCase("/home/pixel/Documents/Cinvestav_2025/Analisis/Chilli/ltpmsi-chilli.imzML", 420.0, 0.1, "Chilli Pepper"),
+    BenchmarkCase("/home/pixel/Documents/Cinvestav_2025/Analisis/imzML_DESI/40TopL,10TopR,30BottomL,20BottomR/40TopL,10TopR,30BottomL,20BottomR-centroid.imzML", 885.5, 0.1, "Colon Cancer Human"),
+    BenchmarkCase("/home/pixel/Documents/Cinvestav_2025/Analisis/imzML_AP_SMALDI/HR2MSImouseurinarybladderS096.imzML", 716.053, 0.1, "Mouse Urinary Bladder"),
+    BenchmarkCase("/home/pixel/Documents/Cinvestav_2025/Analisis/Leafs/CE1_Leaf_R3.imzML", 306.1, 0.1, "Leaf"),
+    BenchmarkCase("/home/pixel/Documents/Cinvestav_2025/Analisis/Liv2_imzML_TIMSConvert-selected/Liv2.imzML", 796.18, 0.1, "Liver Cut"),
+    BenchmarkCase("/home/pixel/Documents/Cinvestav_2025/Analisis/salida_Tims/Stomach_DHB.imzML", 804.3, 0.1, "Mouse Stomach 4GB")
+]
+
+function get_total_file_size_mb(filepath::String)
+    imzml_size = isfile(filepath) ? filesize(filepath) : 0
+    ibd_path = replace(filepath, r"\.(imzML|imzml)$" => ".ibd")
+    ibd_size = isfile(ibd_path) ? filesize(ibd_path) : 0
+    return round((imzml_size + ibd_size) / 1024^2, digits=2)
+end
+
+function flush_memory()
+    GC.gc(true)
+    if Sys.islinux()
+        # Force the OS to reclaim memory from the glibc allocator
+        ccall(:malloc_trim, Int32, (Int32,), 0)
+    end
+end
+
+function run_v3_benchmarks()
+    mkpath(RESULTS_DIR)
+    results = DataFrame()
+
+    println("="^60)
+    println("STARTING ENTERPRISE BENCHMARK SUITE (v3)")
+    println("="^60)
+
+    for case in BENCHMARK_CASES
+        if !isfile(case.filepath)
+            @warn "File not found: $(case.filepath). Skipping."
+            continue
+        end
+
+        println("\n--- Target: $(case.name) ---")
+        file_size = get_total_file_size_mb(case.filepath)
+
+        # ------------------------------------------------------------
+        # 1. JuliaMSI (New Architecture)
+        # ------------------------------------------------------------
+        println("[JuliaMSI - New Engine]")
+        flush_memory()
+        
+        # Load Phase (Metadata Only / Memory Mapping)
+        load_stats_new = @timed OpenMSIData(case.filepath)
+        msi_data = load_stats_new.value
+        load_time_new_s = load_stats_new.time
+        mem_load_new_mb = load_stats_new.bytes / 1024^2
+        
+        # Slicing Phase (High Precision)
+        b_slice_new = @benchmark get_mz_slice($msi_data, $(case.mz_value), $(case.mz_tolerance)) samples=10 seconds=5
+        mean_time_new_s = mean(b_slice_new.times) / 1e9
+        
+        # Calculate Amortized Throughput (10 slices)
+        # Includes the 'Loading Wall' penalty
+        total_time_10_new = load_time_new_s + (10 * mean_time_new_s)
+        amortized_sps_new = 10.0 / total_time_10_new
+        
+        close(msi_data)
+        flush_memory()
+
+        # ------------------------------------------------------------
+        # 2. julia_mzML_imzML (Legacy Architecture)
+        # ------------------------------------------------------------
+        println("[julia_mzML_imzML - Legacy]")
+        
+        load_stats_old = try
+            @timed LoadImzml(case.filepath)
+        catch e
+            @warn "Legacy load failed: $e"
+            (time=Inf, bytes=Inf, value=nothing)
+        end
+        
+        old_data = load_stats_old.value
+        mem_load_old_mb = load_stats_old.bytes / 1024^2
+        
+        load_time_old_s = Inf
+        mean_time_old_s = Inf
+        amortized_sps_old = 0.0
+        
+        if old_data !== nothing
+            b_slice_old = try
+                # Legacy can be extremely slow, constrain it heavily
+                @benchmark GetSlice($old_data, $(case.mz_value), $(case.mz_tolerance)) samples=3 seconds=10
+            catch e
+                @warn "Legacy slice failed: $e"
+                nothing
+            end
+            
+            if b_slice_old !== nothing
+                mean_time_old_s = mean(b_slice_old.times) / 1e9
+                
+                # Accurately reflect the massive load time in the throughput
+                load_time_old_s = load_stats_old.time
+                total_time_10_old = load_time_old_s + (10 * mean_time_old_s)
+                amortized_sps_old = 10.0 / total_time_10_old
+            end
+        end
+        
+        flush_memory()
+
+        # ------------------------------------------------------------
+        # Record
+        # ------------------------------------------------------------
+        push!(results, (
+            Dataset = case.name,
+            FileSize_MB = file_size,
+            
+            # Legacy Metrics
+            Legacy_LoadMem_MB = mem_load_old_mb,
+            Legacy_LoadTime_s = load_stats_old.time,
+            Legacy_Amortized_SPS_10 = amortized_sps_old,
+            
+            # New Metrics
+            New_LoadMem_MB = mem_load_new_mb,
+            New_LoadTime_s = load_time_new_s,
+            New_Amortized_SPS_10 = amortized_sps_new,
+            
+            # Competitive Deltas
+            RAM_Reduction_Pct = isfinite(mem_load_old_mb) ? ((mem_load_old_mb - mem_load_new_mb) / mem_load_old_mb)*100 : NaN,
+            UX_Speedup_Factor = isfinite(amortized_sps_old) ? (amortized_sps_new / amortized_sps_old) : NaN
+        ))
+        
+        println("  > RAM Reduction: $(round(results[end, :RAM_Reduction_Pct], digits=2))%")
+        println("  > Amortized Throughput (10 Slices): $(round(amortized_sps_old, digits=2)) -> $(round(amortized_sps_new, digits=2)) slices/sec")
+    end
+
+    csv_path = joinpath(RESULTS_DIR, "v3_enterprise_benchmarks.csv")
+    CSV.write(csv_path, results)
+    println("\nData saved to $csv_path")
+    
+    plot_v3_results(results)
+    return results
+end
+
+function plot_v3_results(df::DataFrame)
+    if isempty(df) return end
+
+    fig = Figure(size=(1800, 1200), fontsize=24)
+    x_pos = 1:nrow(df)
+    labels = df.Dataset
+
+    # ---------- Plot 1: The RAM Revolution ----------
+    ax1 = Axis(fig[1, 1], 
+               title="Initial RAM Cost (Loading & Mmap)", 
+               ylabel="Memory (MB)",
+               xticks=(x_pos, labels), xticklabelrotation=π/8)
+    
+    barplot!(ax1, x_pos .- 0.2, df.New_LoadMem_MB, color="#10b981", width=0.4, label="JuliaMSI (Mmap Lazy)")
+    barplot!(ax1, x_pos .+ 0.2, df.Legacy_LoadMem_MB, color="#ef4444", width=0.4, label="Legacy (Dense Matrix)")
+    axislegend(ax1, position=:lt)
+
+    # ---------- Plot 2: Throughput Leap ----------
+    ax2 = Axis(fig[1, 2], 
+               title="Amortized Throughput (10 Slices - Higher is Better)", 
+               ylabel="Slices / Second (Inc. Load Time)",
+               xticks=(x_pos, labels), xticklabelrotation=π/8)
+               
+    barplot!(ax2, x_pos .- 0.2, df.New_Amortized_SPS_10, color="#10b981", width=0.4, label="JuliaMSI")
+    barplot!(ax2, x_pos .+ 0.2, df.Legacy_Amortized_SPS_10, color="#ef4444", width=0.4, label="Legacy")
+    axislegend(ax2, position=:lt)
+
+    # ---------- Plot 3: The Scaling Wall (File Size vs Throughput) ----------
+    ax3 = Axis(fig[2, 1:2], 
+               title="Performance Scaling: UX Throughput vs Dataset Size", 
+               xlabel="Dataset File Size (MB)", 
+               ylabel="Amortized Throughput (Slices/sec)")
+               
+    scatterlines!(ax3, df.FileSize_MB, df.New_Amortized_SPS_10, color="#10b981", markersize=15, linewidth=4, label="JuliaMSI Engine")
+    scatterlines!(ax3, df.FileSize_MB, df.Legacy_Amortized_SPS_10, color="#ef4444", markersize=15, linewidth=4, label="Legacy Engine")
+    axislegend(ax3, position=:rt)
+
+    save(joinpath(RESULTS_DIR, "v3_enterprise_dashboard.png"), fig)
+    println("\nDashboards saved successfully.")
+end
+
+# Ensure we process if run as the main script
+if abspath(PROGRAM_FILE) == @__FILE__
+    run_v3_benchmarks()
+end

test/new_benchmark_mmap.jl

@@ -0,0 +1,48 @@
+using BenchmarkTools
+using MSI_src
+using Statistics
+using DataFrames
+
+# ===================================================================
+# HIGH-PRECISION COMPARATIVE SUITE
+# ===================================================================
+
+function run_advanced_benchmark(path, mz, tol)
+    println("\n" * "="^40)
+    println("TARGET: $(basename(path))")
+    println("="^40)
+
+    # 1. NEW LIBRARY: Metadata Load (The "Control Tower" startup)
+    # This measures how fast the Mmap and Cache system works
+    t_load_new = @belapsed OpenMSIData($path)
+    
+    # 2. NEW LIBRARY: Slice Generation (The "Streaming" speed)
+    msi_new = OpenMSIData(path)
+    # We use @benchmark to get a distribution (min, mean, max)
+    b_slice_new = @benchmark get_mz_slice($msi_new, $mz, $tol)
+
+    # --- Metrics Table ---
+    results = DataFrame(
+        Metric = ["Metadata Load", "Slice Gen (Min)", "Slice Gen (Mean)", "Allocations"],
+        JuliaMSI = [
+            "$(round(t_load_new * 1000, digits=2)) ms",
+            "$(round(minimum(b_slice_new.times)/1e6, digits=2)) ms",
+            "$(round(mean(b_slice_new.times)/1e6, digits=2)) ms",
+            "$(b_slice_new.allocs) allocs"
+        ]
+    )
+    
+    println(results)
+    
+    # --- The "Throughput" Test ---
+    # How many slices per second can we handle?
+    throughput_new = 1.0 / mean(b_slice_new.times/1e9)
+    println("\nThroughput: $(round(throughput_new, digits=1)) slices/sec")
+    
+    return results
+end
+
+# Example Run
+@time run_advanced_benchmark("/home/pixel/Documents/Cinvestav_2025/Analisis/imzML_AP_SMALDI/HR2MSImouseurinarybladderS096.imzML", 716.053, 0.1)
+
+# For multithread: julia --threads auto --project=. test/new_benchmark_mmap.jl

test/run_precalculation_example.jl

@@ -16,8 +16,9 @@ using MSI_src
 const TEST_MZML_FILE = ""
 # const TEST_IMZML_FILE = "/home/pixel/Documents/Cinvestav_2025/Analisis/CE4_BF_R1/CE4_BF_R1.imzML"
 #const TEST_IMZML_FILE = "/home/pixel/Documents/Cinvestav_2025/Analisis/salida/Stomach_DHB_uncompressed.imzML"
-const TEST_IMZML_FILE = "/home/pixel/Documents/Cinvestav_2025/Analisis/Thricoderma_etc/Imaging_interaccion_trichoderma_vs_streptomyces.imzML"
+# const TEST_IMZML_FILE = "/home/pixel/Documents/Cinvestav_2025/Analisis/Thricoderma_etc/Imaging_interaccion_trichoderma_vs_streptomyces.imzML"
 #const MASK_ROUTE = "/home/pixel/Documents/Cinvestav_2025/JuliaMSI/public/css/masks/Stomach_DHB_uncompressed.png"
+const TEST_IMZML_FILE = "/home/pixel/Documents/Cinvestav_2025/Analisis/imzML_AP_SMALDI/HR2MSImouseurinarybladderS096.imzML"
 const MASK_ROUTE = ""
 
 #=

test/run_preprocessing.jl

@@ -19,11 +19,11 @@ using MSI_src
 
 # const TEST_FILE = "/home/pixel/Documents/Cinvestav_2025/Analisis/Thricoderma_etc/Imaging_interaccion_trichoderma_vs_streptomyces.imzML"
 # const TEST_FILE = "/home/pixel/Documents/Cinvestav_2025/Analisis/set de datos MS/Atropina_tuneo_fraq_20ev.mzML"
-const TEST_FILE = "/home/pixel/Documents/Cinvestav_2025/Analisis/salida/Stomach_DHB_uncompressed.imzML"
+# const TEST_FILE = "/home/pixel/Documents/Cinvestav_2025/Analisis/salida/Stomach_DHB_uncompressed.imzML"
-# const TEST_FILE = "/home/pixel/Documents/Cinvestav_2025/Analisis/imzML_AP_SMALDI/HR2MSImouseurinarybladderS096.imzML"
+const TEST_FILE = "/home/pixel/Documents/Cinvestav_2025/Analisis/imzML_AP_SMALDI/HR2MSImouseurinarybladderS096.imzML"
 
-const MASK_ROUTE = "/home/pixel/Documents/Cinvestav_2025/JuliaMSI/public/css/masks/Stomach_DHB_uncompressed.png"
+# const MASK_ROUTE = "/home/pixel/Documents/Cinvestav_2025/JuliaMSI/public/css/masks/Stomach_DHB_uncompressed.png"
-# const MASK_ROUTE = "" 
+const MASK_ROUTE = "" 
 
 const OUTPUT_DIR = "./test/results/preprocessing_results"
 

test/run_tests.jl

@@ -31,9 +31,10 @@ using MSI_src
 # --- Test Case 1: Standard .mzML file ---
 # A regular, non-imaging mzML file.
 # const TEST_MZML_FILE = "/home/pixel/Documents/Cinvestav_2025/Analisis/mzML/T9_A1.mzML"
-const TEST_MZML_FILE = "/home/pixel/Documents/Cinvestav_2025/Analisis/CE4_BF_R1/CE4_BF_R1.mzML"
+# const TEST_MZML_FILE = "/home/pixel/Documents/Cinvestav_2025/Analisis/CE4_BF_R1/CE4_BF_R1.mzML"
 # const TEST_MZML_FILE = "/home/pixel/Documents/Cinvestav_2025/Analisis/Imaging_paper_spray/Imaging_paper_spray.mzML"
 # const TEST_MZML_FILE = "/home/pixel/Documents/Cinvestav_2025/Analisis/Imaging prueba Roya 1/Roya.mzML"
+const TEST_MZML_FILE = "/home/pixel/Documents/Cinvestav_2025/Analisis/mzML"
 const SPECTRUM_TO_PLOT = 1 # Which spectrum to plot from the file
 
 # --- Test Case 2: .mzML + Sync File for Conversion ---
@@ -58,16 +59,17 @@ const CONVERSION_TARGET_IMZML = "test/results/converted_mzml.imzML"
 # const TEST_IMZML_FILE = "/home/pixel/Documents/Cinvestav_2025/Analisis/Imaging prueba Roya 1/royaimg.imzML"
 # const TEST_IMZML_FILE = "/home/pixel/Documents/Cinvestav_2025/Analisis/salida/ltpmsi-chilli.imzML" # centroid aparently?
 # const TEST_IMZML_FILE = "/home/pixel/Documents/Cinvestav_2025/Analisis/salida/Stomach_DHB_compressed.imzML" # centroid compressed
-const TEST_IMZML_FILE = "/home/pixel/Documents/Cinvestav_2025/Analisis/salida/Stomach_DHB_uncompressed.imzML" # centroid
+# const TEST_IMZML_FILE = "/home/pixel/Documents/Cinvestav_2025/Analisis/salida/Stomach_DHB_uncompressed.imzML" # centroid
+const TEST_IMZML_FILE = "/home/pixel/Documents/Cinvestav_2025/Analisis/imzML_AP_SMALDI/HR2MSImouseurinarybladderS096.imzML"
 # The m/z value to use for creating an image slice.
 # const MZ_VALUE_FOR_SLICE = 309.06 # BF
-# const MZ_VALUE_FOR_SLICE = 896.0 # HR2MSI
+const MZ_VALUE_FOR_SLICE = 896.0 # HR2MSI
 # const MZ_VALUE_FOR_SLICE = 76.03 # I PS
 # const MZ_VALUE_FOR_SLICE = 313 # ROYA
-const MZ_VALUE_FOR_SLICE = 100 # advanced processing
+# const MZ_VALUE_FOR_SLICE = 100 # advanced processing
 # const MZ_TOLERANCE = 0.1
 # const MZ_TOLERANCE = 1
-const MZ_TOLERANCE = 0.1
+const MZ_TOLERANCE = 0.2
 
 # Coordinates to plot a specific spectrum from imzML
 const COORDS_TO_PLOT = (50, 50) # Example coordinates (X, Y)

test/runtests.jl

@@ -0,0 +1,58 @@
+using Test
+
+# We need to simulate loading the MSI_src module directly
+include("../src/MSI_src.jl")
+using .MSI_src
+
+@testset "JuliaMSI Core Systems" begin
+    @testset "Basic Instantiation" begin
+        # 1. Test basic metadata structure initialization
+        meta = SpectrumMetadata(
+            Int32(1), Int32(2), 
+            "test_id", 
+            :sample, 
+            CENTROID, 
+            SpectrumAsset(Float64, false, Int64(0), 100, :mz, 0.0, 0.0), 
+            SpectrumAsset(Float32, true, Int64(100), 50, :intensity, 0.0, 0.0)
+        )
+        @test meta.x == 1
+        @test meta.y == 2
+        @test meta.mode == CENTROID
+        @test meta.mz_asset.format == Float64
+        @test meta.int_asset.format == Float32
+        @test meta.int_asset.is_compressed == true
+
+        # 2. Test cache pool and MSIData structural stability
+        source = MzMLSource([], Float64, Float32, nothing) # Empty source for structural testing
+        
+        # Test constructor doesn't throw
+        # Constructor signature: (source, metadata, instrument_meta, dims, coordinate_map, cache_size)
+        msi_data = MSIData(
+            source, 
+            [meta], 
+            nothing, # instrument_meta
+            (100, 100), # dims
+            nothing, # coord map
+            10 # cache size
+        )
+        
+        @test msi_data.image_dims == (100, 100)
+        @test length(msi_data.spectra_metadata) == 1
+        @test msi_data.cache_size == 10
+    end
+
+    @testset "Buffer & Cache Subsystems" begin
+        pool = SimpleBufferPool()
+        # Test basic retrieval
+        buf = get_buffer!(pool, 1024)
+        @test length(buf) == 1024
+        
+        # Test release
+        release_buffer!(pool, buf)
+        @test length(pool.buffers[1024]) == 1
+        
+        # Test reuse
+        buf2 = get_buffer!(pool, 1024)
+        @test buf === buf2 # Should return the EXACT same buffer object
+    end
+end

test/test_buffer_pool_concurrency.jl

@@ -0,0 +1,32 @@
+using MSI_src
+using Test
+using Base.Threads
+
+@testset "SimpleBufferPool Concurrency Stress Test" begin
+    pool = MSI_src.SimpleBufferPool()
+    n_threads = Threads.nthreads()
+    n_iterations = 1000
+    buffer_size = 1024
+
+    println("Running buffer pool stress test with $n_threads threads...")
+
+    # Parallel stress test
+    Threads.@threads for i in 1:(n_threads * n_iterations)
+        # get_buffer! and release_buffer! are now thread-safe
+        buf = MSI_src.get_buffer!(pool, buffer_size)
+        
+        # Simulate some work
+        fill!(buf, UInt8(i % 256))
+        
+        MSI_src.release_buffer!(pool, buf)
+    end
+
+    # After stress test, the dictionary should be coherent
+    sizes = collect(keys(pool.buffers))
+    if !isempty(sizes)
+        @test buffer_size ∈ sizes
+        @test length(pool.buffers[buffer_size]) <= pool.max_pool_size
+    end
+    
+    println("Buffer pool stress test PASSED.")
+end

test/test_pipeline_concurrency.jl

@@ -0,0 +1,146 @@
+
+using Pkg
+Pkg.activate(joinpath(@__DIR__, ".."))
+using MSI_src
+using Test
+using Base.Threads
+
+# --- Test Configuration ---
+const TEST_FILE = "/home/pixel/Documents/Cinvestav_2025/Analisis/Stomach/Stomach_DHB_uncompressed.imzML"
+
+println("Starting concurrency test with file: $TEST_FILE")
+
+if !isfile(TEST_FILE)
+    error("Test file not found: $TEST_FILE")
+end
+
+# --- Simulation of the Concurrency Issue ---
+
+# Global "UI State" variable simulates the app's msi_data
+global_msi_data = nothing
+const GLOBAL_LOCK = ReentrantLock()
+
+function ui_load_file(path)
+    global global_msi_data
+    lock(GLOBAL_LOCK) do
+        if global_msi_data !== nothing
+            close(global_msi_data)
+        end
+        println("UI: Loading file...")
+        global_msi_data = OpenMSIData(path)
+        println("UI: File loaded.")
+    end
+end
+
+function ui_close_file()
+    global global_msi_data
+    lock(GLOBAL_LOCK) do
+        if global_msi_data !== nothing
+            println("UI: Closing file...")
+            close(global_msi_data)
+            global_msi_data = nothing
+            println("UI: File closed.")
+        end
+    end
+end
+
+# Mock pipeline function that mimics `run_full_pipeline` in app.jl
+# Crucially, it uses its own local `pipeline_msi_data` as per the fix
+function run_mock_pipeline_isolated(path)
+    println("Pipeline: Starting isolated pipeline...")
+    
+    # 1. Open ISOLATED instance
+    pipeline_msi_data = OpenMSIData(path)
+    println("Pipeline: Opened isolated MSIData instance.")
+    
+    try
+        # 2. Simulate reading spectra
+        indices = 1:min(100, length(pipeline_msi_data.spectra_metadata)) # Read first 100 spectra
+        
+        # Artificial delay to allow "user" interaction
+        sleep(0.5) 
+        
+        println("Pipeline: Reading spectra...")
+        Threads.@threads for i in indices
+            # Use local instance
+            mz, int = GetSpectrum(pipeline_msi_data, i)
+            # Simulate processing work
+            sum(int) 
+        end
+        println("Pipeline: Finished reading spectra successfully.")
+        return true
+    catch e
+        println("Pipeline: CRASHED with error: $e")
+        return false
+    finally
+        close(pipeline_msi_data)
+        println("Pipeline: Closed isolated MSIData instance.")
+    end
+end
+
+# Mock pipeline that uses GLOBAL instance (The BUGGY version)
+function run_mock_pipeline_buggy()
+    println("Buggy Pipeline: Starting...")
+    # Uses global_msi_data directly
+    
+    try
+        global global_msi_data
+        if global_msi_data === nothing
+            println("Buggy Pipeline: No data loaded!")
+            return false
+        end
+        
+        local_ref = global_msi_data # Still points to same object
+        
+        indices = 1:min(100, length(local_ref.spectra_metadata))
+        
+        sleep(0.5)
+        
+        println("Buggy Pipeline: Reading spectra from shared object...")
+        Threads.@threads for i in indices
+            # This will fail if ui_close_file() happens concurrently
+            mz, int = GetSpectrum(local_ref, i)
+            sum(int)
+        end
+        println("Buggy Pipeline: Success (Unexpected if concurrency worked)")
+        return true
+    catch e
+        println("Buggy Pipeline: CRASHED as expected: $e")
+        return false
+    end
+end
+
+
+# --- execute checks ---
+
+@testset "Concurrency Crash Fix Verification" begin
+    
+    # 1. Setup: Load file initially
+    ui_load_file(TEST_FILE)
+    
+    # 2. Test the FIX: Isolated Pipeline
+    println("\n--- Testing Fixed (Isolated) Pipeline ---")
+    
+    t_pipeline = Threads.@spawn run_mock_pipeline_isolated(TEST_FILE)
+    
+    # Simulate user closing/reloading file while pipeline runs
+    sleep(0.2) 
+    ui_close_file()
+    
+    # Wait for pipeline
+    success = fetch(t_pipeline)
+    @test success == true
+    println("Fixed pipeline result: ", success ? "PASSED" : "FAILED")
+
+
+    # 3. Test the BUG: Global Pipeline (Optional, to prove it crashes without fix)
+    # Uncomment to verify the bug exists if needed, but we assume it does based on user report.
+    # println("\n--- Testing Buggy (Shared) Pipeline ---")
+    # ui_load_file(TEST_FILE)
+    # t_buggy = Threads.@spawn run_mock_pipeline_buggy()
+    # sleep(0.2)
+    # ui_close_file()
+    # buggy_success = fetch(t_buggy)
+    # println("Buggy pipeline result: ", buggy_success ? "PASSED (No crash?)" : "FAILED (Crashed as expected)")
+    
+end

test/test_streaming_pipeline.jl

@@ -0,0 +1,124 @@
+#!/usr/bin/env julia
+# test/test_streaming_pipeline.jl
+# ============================================================================
+# Validation test for Sprint 2: The Streaming Pipeline
+#
+# This test exercises process_dataset! against the HR2MSI mouse bladder
+# dataset and verifies:
+#   1. Correct sparse matrix creation
+#   2. Non-zero peak population
+#   3. RAM savings vs dense equivalent
+#   4. Allocation count and throughput
+# ============================================================================
+
+using Pkg
+Pkg.activate(".")
+
+using MSI_src
+using SparseArrays
+
+# =============================================================================
+# Configuration
+# =============================================================================
+
+const IMZML_PATH = "/home/pixel/Documents/Cinvestav_2025/Analisis/imzML_AP_SMALDI/HR2MSImouseurinarybladderS096.imzML"
+
+function main()
+    println("=" ^ 60)
+    println("SPRINT 2: Streaming Pipeline Validation")
+    println("=" ^ 60)
+    
+    if !isfile(IMZML_PATH)
+        println("SKIPPED: Dataset not found at $IMZML_PATH")
+        return
+    end
+    
+    # --- 1. Load dataset ---
+    println("\n--- Step 1: Loading dataset ---")
+    data = OpenMSIData(IMZML_PATH)
+    println("Loaded: $(length(data.spectra_metadata)) spectra")
+    
+    # --- 2. Configure the streaming pipeline ---
+    println("\n--- Step 2: Configuring pipeline ---")
+    config = PipelineConfig(
+        steps = [
+            StreamingStep(:baseline_correction, Dict{Symbol,Any}(:method => :snip, :iterations => 50)),
+            StreamingStep(:normalization, Dict{Symbol,Any}(:method => :tic)),
+            StreamingStep(:peak_picking, Dict{Symbol,Any}(
+                :method => :profile,
+                :snr_threshold => 3.0,
+                :half_window => 10,
+                :min_peak_prominence => 0.1,
+                :merge_peaks_tolerance => 0.002
+            )),
+        ],
+        num_bins = 2000,
+        frequency_threshold = 0.01  # Bins must appear in at least 1% of spectra
+    )
+    println("Steps: $(join([s.name for s in config.steps], " → "))")
+    println("Bins: $(config.num_bins), Frequency threshold: $(config.frequency_threshold)")
+    
+    # --- 3. Run the streaming pipeline ---
+    println("\n--- Step 3: Running streaming pipeline ---")
+    stats = @timed begin
+        feature_matrix, bin_centers = process_dataset!(data, config)
+    end
+    
+    feature_matrix = stats.value[1]
+    bin_centers = stats.value[2]
+    
+    println("\n--- Results ---")
+    println("  Feature matrix size: $(size(feature_matrix))")
+    println("  Non-zeros: $(nnz(feature_matrix))")
+    println("  Bin centers: $(length(bin_centers))")
+    println("  Time: $(round(stats.time, digits=2))s")
+    println("  Allocations: $(stats.bytes ÷ 1_000_000) MB")
+    println("  GC time: $(round(stats.gctime, digits=2))s")
+    
+    # --- 4. Validate ---
+    println("\n--- Step 4: Validation ---")
+    
+    passed = true
+    
+    # Check matrix dimensions
+    if size(feature_matrix, 1) > 0 && size(feature_matrix, 2) > 0
+        println("  ✓ Matrix has valid dimensions")
+    else
+        println("  ✗ Matrix has invalid dimensions: $(size(feature_matrix))")
+        passed = false
+    end
+    
+    # Check non-zeros
+    if nnz(feature_matrix) > 0
+        println("  ✓ Matrix has $(nnz(feature_matrix)) non-zero entries")
+    else
+        println("  ✗ Matrix is completely empty")
+        passed = false
+    end
+    
+    # Check sparsity savings
+    dense_mb = size(feature_matrix, 1) * size(feature_matrix, 2) * 8 / 1e6
+    sparse_mb = nnz(feature_matrix) * 16 / 1e6  # index + value per entry
+    if dense_mb > 0
+        savings = (1.0 - sparse_mb / dense_mb) * 100
+        println("  ✓ RAM savings: $(round(savings, digits=1))% ($(round(sparse_mb, digits=1)) MB vs $(round(dense_mb, digits=1)) MB dense)")
+    end
+    
+    # Check bin centers alignment
+    if length(bin_centers) == size(feature_matrix, 1)
+        println("  ✓ Bin centers match matrix rows")
+    else
+        println("  ✗ Bin center count ($(length(bin_centers))) != matrix rows ($(size(feature_matrix, 1)))")
+        passed = false
+    end
+    
+    println("\n" * "=" ^ 60)
+    if passed
+        println("ALL VALIDATIONS PASSED ✓")
+    else
+        println("SOME VALIDATIONS FAILED ✗")
+    end
+    println("=" ^ 60)
+end
+
+@time main()