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

@@ -8,3 +8,5 @@ R_original_scripts/
 test/results/
 test/binarization_results/
 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"]

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

@@ -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.

app.jl

@@ -29,6 +29,9 @@ if !@isdefined(increment_image)
     include("./julia_imzML_visual.jl")
 end
 
+const global_msi_data = Ref{Union{MSIData, Nothing}}(nothing)
+
+
 # --- Memory Validation Logging ---
 if get(ENV, "GENIE_ENV", "dev") != "prod"
     function get_rss_mb()
@@ -60,7 +63,7 @@ if get(ENV, "GENIE_ENV", "dev") != "prod"
         println("--- MEMORY LOG [$(context)] ---")
         println("  Timestamp: $(now())")
         println("  Process RSS: $(rss_mb) MB")
-        println("  msi_data size: $(msi_data_size_mb) MB")
+        println("  global_msi_data[] size: $(msi_data_size_mb) MB")
         println("  Cumulative GC time: $(gc_time_s) s")
         println("--------------------------")
     end
@@ -193,6 +196,15 @@ macro ui_log(message, level="INFO", log_entries)
 end
 =#
 
+# --- CRITICAL: Disable Stipple's session-to-disk persistence ---
+# Stipple's ModelStorage registers on(field) handlers that serialize the ENTIRE
+# ReactiveModel to disk via GenieSessionFileSession on every UI state change.
+# With 253 reactive variables including multi-MB Plotly traces, this generates
+# gigabytes of orphaned session files in /tmp/jl_XXXXXX, exhausting disk space.
+# For a single-user desktop application, session persistence is unnecessary.
+Stipple.enable_model_storage(false)
+Core.eval(Stipple, :(sesstoken() = "")) # Prevent ErrorException("Model storage is disabled") during layout render
+
 # Reactive code to make the UI interactive
 @app begin
     # == Notification & Logs ==
@@ -476,7 +488,7 @@ end
 
     # == DATA MANAGEMENT VARIABLES ==
     # Centralized MSIData object
-    @out msi_data::Union{MSIData, Nothing} = nothing
+    # global_msi_data[] is now global to avoid Genie Session memory leaks
 
     # Image file management
     @out text_nmass="" # For specific mass charge image creation
@@ -661,7 +673,7 @@ end
         
         try
             # 1. Clear large data objects explicitly
-            msi_data = nothing
+            global_msi_data[] = nothing
             feature_matrix_result = nothing
             bin_info_result = nothing
             
@@ -717,29 +729,41 @@ end
     sure the file can be processed by later steps like mainProcess 
     =#
     @onbutton btnSearch begin
-        is_processing = true
-        push!(__model__)
+        # 0. Robustness Guard: Prevent double-trigger during processing
+        if is_processing
+            println("DEBUG: btnSearch ignored because another process is already running.")
+            return
+        end
+        btnSearch = false # Manual reset of the trigger
+
+        # 1. Grab the file path from the main task
         picked_route = pick_file(; filterlist="imzML,imzml,mzML,mzml")
         
         if isnothing(picked_route) || isempty(picked_route)
-            is_processing = false
             return
         end
 
+        # 2. Update reactive state synchronously
+        is_processing = true
+        msg = "Opening file: $(basename(picked_route))..."
+        SpectraEnabled = false
+        btnMetadataDisable = true
+        push!(__model__)
+
+        # 3. Spawn background computational thread
+        Threads.@spawn begin
+            try
                 # --- Close previous dataset if one is open ---
-        if msi_data !== nothing
-            println("DEBUG: Closing previously loaded dataset before opening new one: $(basename(full_route))")
-            close(msi_data)
-            msi_data = nothing
+                if global_msi_data[] !== nothing
+                    println("DEBUG: Closing previously loaded dataset before opening new one...")
+                    close(global_msi_data[])
+                    global_msi_data[] = nothing
                     GC.gc()
                     if Sys.islinux()
                         ccall(:malloc_trim, Int32, (Int32,), 0)
                     end
                 end
 
-        msg = "Opening file: $(basename(picked_route))..."
-
-        try
                 dataset_name = replace(basename(picked_route), r"(\.(imzML|imzml|mzML|mzml))$"i => "")
             registry = load_registry(registry_path)
             existing_entry = get(registry, dataset_name, nothing)
@@ -757,8 +781,8 @@ end
                 dims = parse.(Int, split(dims_str, " x "))
                 imgWidth, imgHeight = dims[1], dims[2]
 
-                msi_data = nothing # Ensure data is not held in memory
-                log_memory_usage("Fast Load (msi_data cleared)", msi_data)
+                global_msi_data[] = nothing # Ensure data is not held in memory
+                log_memory_usage("Fast Load (global_msi_data[] cleared)", global_msi_data[])
                 btnMetadataDisable = false
                 SpectraEnabled = true
                 selected_folder_main = dataset_name
@@ -1000,10 +1024,10 @@ end
             image_available_folders = deepcopy(img_folders)
 
             selected_folder_main = dataset_name
-            msi_data = loaded_data
+            global_msi_data[] = loaded_data
 
             # Determine plot mode from loaded data
-            df = msi_data.spectrum_stats_df
+            df = global_msi_data[].spectrum_stats_df
             if df !== nothing && "Mode" in names(df)
                 profile_count = count(==(MSI_src.PROFILE), df.Mode)
                 total_count = length(df.Mode)
@@ -1013,7 +1037,7 @@ end
                 last_plot_mode = "lines" # Default
             end
 
-            log_memory_usage("Full Load", msi_data)
+            log_memory_usage("Full Load", global_msi_data[])
 
             eTime = round(time() - sTime, digits=3)
             msg = "Active file loaded in $(eTime) seconds. Dataset '$(dataset_name)' is ready for analysis."
@@ -1021,18 +1045,21 @@ end
             SpectraEnabled = true
             
             catch e
-            msi_data = nothing
+                global_msi_data[] = nothing
                 msg = "Error loading active file: $e"
                 warning_msg = true
                 SpectraEnabled = false
                 btnMetadataDisable = true
                 @error "File loading failed" exception=(e, catch_backtrace())
+                push!(__model__) # Force sending error back to UI immediately
             finally
                 GC.gc() # Trigger garbage collection
                 if Sys.islinux()
                     ccall(:malloc_trim, Int32, (Int32,), 0) # Ensure Julia returns the freed memory to OS
                 end
                 is_processing = false
+                push!(__model__) # Clear spinner loop 
+            end
         end
     end
 
@@ -1232,6 +1259,11 @@ end
     This reactive handler job is to run the full preprocessing pipeline on the selected dataset.
     =#
     @onbutton run_full_pipeline begin
+        if is_processing
+            println("DEBUG: run_full_pipeline ignored because another process is already running.")
+            return
+        end
+        run_full_pipeline = false # Manual reset
         is_processing = true
         push!(__model__)
         overall_progress = 0.0
@@ -1259,9 +1291,9 @@ end
             end
             target_path = entry["source_path"]
 
-            # Ensure msi_data is for the currently selected file and load if needed
+            # Ensure global_msi_data[] is for the currently selected file and load if needed
             # NOTE: For the pipeline, we will open a DEDICATED instance to avoid race conditions
-            # with the global msi_data used for plotting/interactive exploration.
+            # with the global global_msi_data[] used for plotting/interactive exploration.
             println("DEBUG: Opening isolated MSIData instance for pipeline stability...")
             pipeline_msi_data = OpenMSIData(target_path)
 
@@ -1434,7 +1466,11 @@ end
             # --- 2. Parameter Assembly with Validation ---
             current_pipeline_step = "Configuring parameters..."
             println("DEBUG: Configuring parameters and validating enabled steps...")
-            ref_peaks = Dict{Float64, String}(p["mz"] => p["label"] for p in reference_peaks_list)
+            ref_peaks = Dict{Float64, String}(
+                parse(Float64, string(p["mz"])) => string(get(p, "label", ""))
+                for p in reference_peaks_list 
+                if tryparse(Float64, string(p["mz"])) !== nothing
+            )
 
             final_params = Dict{Symbol, Dict{Symbol, Any}}()
             validation_errors = String[]
@@ -1682,9 +1718,9 @@ end
                     
                     # Determine plot mode for this specific spectrum
                     spectrum_mode_for_plot = "lines" # Default to lines
-                    if msi_data.spectrum_stats_df !== nothing && "Mode" in names(msi_data.spectrum_stats_df)
-                        if selected_spectrum_id_for_plot > 0 && selected_spectrum_id_for_plot <= length(msi_data.spectrum_stats_df.Mode)
-                            mode = msi_data.spectrum_stats_df.Mode[selected_spectrum_id_for_plot]
+                    if global_msi_data[].spectrum_stats_df !== nothing && "Mode" in names(global_msi_data[].spectrum_stats_df)
+                        if selected_spectrum_id_for_plot > 0 && selected_spectrum_id_for_plot <= length(global_msi_data[].spectrum_stats_df.Mode)
+                            mode = global_msi_data[].spectrum_stats_df.Mode[selected_spectrum_id_for_plot]
                             if mode == MSI_src.CENTROID
                                 spectrum_mode_for_plot = "stem"
                             end
@@ -1935,7 +1971,7 @@ end
     @onbutton recalculate_suggestions_btn begin
         is_processing = true
         push!(__model__)
-        if msi_data === nothing
+        if global_msi_data[] === nothing
             msg = "Please load a file first."
             warning_msg = true
             return
@@ -1944,8 +1980,12 @@ end
         try
             msg = "Recalculating suggestions..."
             
-            ref_peaks = Dict(p["mz"] => p["label"] for p in reference_peaks_list)
-            recommended_params = main_precalculation(msi_data, reference_peaks=ref_peaks)
+            ref_peaks = Dict{Float64, String}(
+                parse(Float64, string(p["mz"])) => string(get(p, "label", ""))
+                for p in reference_peaks_list 
+                if tryparse(Float64, string(p["mz"])) !== nothing
+            )
+            recommended_params = main_precalculation(global_msi_data[], reference_peaks=ref_peaks)
             
 
             for (step_name, params) in recommended_params
@@ -2340,6 +2380,11 @@ end
     end
     
     @onbutton mainProcess @time begin
+        if is_processing
+            println("DEBUG: mainProcess ignored because another process is already running.")
+            return
+        end
+        mainProcess = false # Manual reset
         # --- UI State Update ---
         overall_progress = 0.0
         progress_message = "Preparing batch process..."
@@ -2414,9 +2459,23 @@ end
                 else
                     push!(newly_created_folders, replace(basename(file_path), r"\.imzML$"i => ""))
                 end
+                
+                # --- Periodic Memory Reclamation ---
+                # Force GC and return memory to the OS after processing each large file
+                GC.gc(true)
+                if Sys.islinux()
+                    ccall(:malloc_trim, Cint, (Cint,), 0)
+                end
+                
                 current_step += 1
             end
             
+            # Final memory cleanup for the batch
+            GC.gc(true)
+            if Sys.islinux()
+                ccall(:malloc_trim, Cint, (Cint,), 0)
+            end
+
             # --- 3. Final Report ---
             total_time_end = round(time() - total_time_start, digits=3)
             
@@ -2550,21 +2609,21 @@ end
                 return
             end
 
-            if msi_data === nothing || full_route != target_path
-                if msi_data !== nothing
-                    close(msi_data)
+            if global_msi_data[] === nothing || full_route != target_path
+                if global_msi_data[] !== nothing
+                    close(global_msi_data[])
                 end
                 msg = "Reloading $(basename(target_path)) for analysis..."
                 full_route = target_path
-                msi_data = OpenMSIData(target_path)
+                global_msi_data[] = OpenMSIData(target_path)
                 if haskey(get(entry, "metadata", Dict()), "global_min_mz") && entry["metadata"]["global_min_mz"] !== nothing
                     raw_min = entry["metadata"]["global_min_mz"]
                     raw_max = entry["metadata"]["global_max_mz"]
                     min_val = isa(raw_min, Dict) ? get(raw_min, "value", raw_min) : raw_min
                     max_val = isa(raw_max, Dict) ? get(raw_max, "value", raw_max) : raw_max
-                    set_global_mz_range!(msi_data, convert(Float64, min_val), convert(Float64, max_val))
+                    set_global_mz_range!(global_msi_data[], convert(Float64, min_val), convert(Float64, max_val))
                 else
-                    precompute_analytics(msi_data)
+                    precompute_analytics(global_msi_data[])
                 end
             end
 
@@ -2577,7 +2636,7 @@ end
                 end
             end
 
-            plotdata, plotlayout, xSpectraMz, ySpectraMz = meanSpectrumPlot(msi_data, selected_folder_main, mask_path=mask_path_for_plot)
+            plotdata, plotlayout, xSpectraMz, ySpectraMz = meanSpectrumPlot(global_msi_data[], selected_folder_main, mask_path=mask_path_for_plot)
             plotdata_before = plotdata
             plotlayout_before = plotlayout
             last_plot_type = "mean"
@@ -2585,7 +2644,7 @@ end
             fTime = time()
             eTime = round(fTime - sTime, digits=3)
             msg = "Plot loaded in $(eTime) seconds"
-            log_memory_usage("Mean Plot Generated", msi_data)
+            log_memory_usage("Mean Plot Generated", global_msi_data[])
         catch e
             msg = "Could not generate mean spectrum plot: $e"
             warning_msg = true
@@ -2639,21 +2698,21 @@ end
                 return
             end
 
-            if msi_data === nothing || full_route != target_path
-                if msi_data !== nothing
-                    close(msi_data)
+            if global_msi_data[] === nothing || full_route != target_path
+                if global_msi_data[] !== nothing
+                    close(global_msi_data[])
                 end
                 msg = "Reloading $(basename(target_path)) for analysis..."
                 full_route = target_path
-                msi_data = OpenMSIData(target_path)
+                global_msi_data[] = OpenMSIData(target_path)
                 if haskey(get(entry, "metadata", Dict()), "global_min_mz") && entry["metadata"]["global_min_mz"] !== nothing
                     raw_min = entry["metadata"]["global_min_mz"]
                     raw_max = entry["metadata"]["global_max_mz"]
                     min_val = isa(raw_min, Dict) ? get(raw_min, "value", raw_min) : raw_min
                     max_val = isa(raw_max, Dict) ? get(raw_max, "value", raw_max) : raw_max
-                    set_global_mz_range!(msi_data, convert(Float64, min_val), convert(Float64, max_val))
+                    set_global_mz_range!(global_msi_data[], convert(Float64, min_val), convert(Float64, max_val))
                 else
-                    precompute_analytics(msi_data)
+                    precompute_analytics(global_msi_data[])
                 end
             end
             local mask_path_for_plot::Union{String, Nothing} = nothing
@@ -2665,7 +2724,7 @@ end
                 end
             end
 
-            plotdata, plotlayout, xSpectraMz, ySpectraMz = sumSpectrumPlot(msi_data, selected_folder_main, mask_path=mask_path_for_plot)
+            plotdata, plotlayout, xSpectraMz, ySpectraMz = sumSpectrumPlot(global_msi_data[], selected_folder_main, mask_path=mask_path_for_plot)
             plotdata_before = plotdata
             plotlayout_before = plotlayout
             last_plot_type = "sum"
@@ -2673,7 +2732,7 @@ end
             fTime = time()
             eTime = round(fTime - sTime, digits=3)
             msg = "Total plot loaded in $(eTime) seconds"
-            log_memory_usage("Sum Plot Generated", msi_data)
+            log_memory_usage("Sum Plot Generated", global_msi_data[])
         catch e
             msg = "Could not generate total spectrum plot: $e"
             warning_msg = true
@@ -2730,21 +2789,21 @@ end
                 return
             end
 
-            if msi_data === nothing || full_route != target_path
-                if msi_data !== nothing
-                    close(msi_data)
+            if global_msi_data[] === nothing || full_route != target_path
+                if global_msi_data[] !== nothing
+                    close(global_msi_data[])
                 end
                 msg = "Reloading $(basename(target_path)) for analysis..."
                 full_route = target_path
-                msi_data = OpenMSIData(target_path)
+                global_msi_data[] = OpenMSIData(target_path)
                 if haskey(get(entry, "metadata", Dict()), "global_min_mz") && entry["metadata"]["global_min_mz"] !== nothing
                     raw_min = entry["metadata"]["global_min_mz"]
                     raw_max = entry["metadata"]["global_max_mz"]
                     min_val = isa(raw_min, Dict) ? get(raw_min, "value", raw_min) : raw_min
                     max_val = isa(raw_max, Dict) ? get(raw_max, "value", raw_max) : raw_max
-                    set_global_mz_range!(msi_data, convert(Float64, min_val), convert(Float64, max_val))
+                    set_global_mz_range!(global_msi_data[], convert(Float64, min_val), convert(Float64, max_val))
                 else
-                    precompute_analytics(msi_data)
+                    precompute_analytics(global_msi_data[])
                 end
             end
 
@@ -2759,7 +2818,7 @@ end
 
             # Convert to positive coordinates for processing
             y_positive = yCoord < 0 ? abs(yCoord) : yCoord
-            plotdata, plotlayout, xSpectraMz, ySpectraMz, spectrum_id = xySpectrumPlot(msi_data, xCoord, y_positive, imgWidth, imgHeight, selected_folder_main, mask_path=mask_path_for_plot)
+            plotdata, plotlayout, xSpectraMz, ySpectraMz, spectrum_id = xySpectrumPlot(global_msi_data[], xCoord, y_positive, imgWidth, imgHeight, selected_folder_main, mask_path=mask_path_for_plot)
             plotdata_before = plotdata
             plotlayout_before = plotlayout
             last_plot_type = "single"
@@ -2800,7 +2859,7 @@ end
             fTime = time()
             eTime = round(fTime - sTime, digits=3)
             msg = "Plot loaded in $(eTime) seconds"
-            log_memory_usage("XY Plot Generated", msi_data)
+            log_memory_usage("XY Plot Generated", global_msi_data[])
         catch e
             msg = "Could not retrieve spectrum: $e"
             warning_msg = true
@@ -2845,21 +2904,21 @@ end
                 return
             end
 
-            if msi_data === nothing || full_route != target_path
-                if msi_data !== nothing
-                    close(msi_data)
+            if global_msi_data[] === nothing || full_route != target_path
+                if global_msi_data[] !== nothing
+                    close(global_msi_data[])
                 end
                 msg = "Reloading $(basename(target_path)) for analysis..."
                 full_route = target_path
-                msi_data = OpenMSIData(target_path)
+                global_msi_data[] = OpenMSIData(target_path)
                 if haskey(get(entry, "metadata", Dict()), "global_min_mz") && entry["metadata"]["global_min_mz"] !== nothing
                     raw_min = entry["metadata"]["global_min_mz"]
                     raw_max = entry["metadata"]["global_max_mz"]
                     min_val = isa(raw_min, Dict) ? get(raw_min, "value", raw_min) : raw_min
                     max_val = isa(raw_max, Dict) ? get(raw_max, "value", raw_max) : raw_max
-                    set_global_mz_range!(msi_data, convert(Float64, min_val), convert(Float64, max_val))
+                    set_global_mz_range!(global_msi_data[], convert(Float64, min_val), convert(Float64, max_val))
                 else
-                    precompute_analytics(msi_data)
+                    precompute_analytics(global_msi_data[])
                 end
             end
 
@@ -2873,7 +2932,7 @@ end
             end
 
             # Call the new nSpectrumPlot function
-            plotdata, plotlayout, xSpectraMz, ySpectraMz, spectrum_id = nSpectrumPlot(msi_data, idSpectrum, selected_folder_main, mask_path=mask_path_for_plot)
+            plotdata, plotlayout, xSpectraMz, ySpectraMz, spectrum_id = nSpectrumPlot(global_msi_data[], idSpectrum, selected_folder_main, mask_path=mask_path_for_plot)
             plotdata_before = plotdata
             plotlayout_before = plotlayout
             last_plot_type = "single"
@@ -2883,7 +2942,7 @@ end
             fTime = time()
             eTime = round(fTime - sTime, digits=3)
             msg = "Plot loaded in $(eTime) seconds"
-            log_memory_usage("nSpectrum Plot Generated", msi_data)
+            log_memory_usage("nSpectrum Plot Generated", global_msi_data[])
         catch e
             msg = "Could not retrieve spectrum: $e"
             warning_msg = true
@@ -3164,7 +3223,7 @@ end
 
     # This handler will now correctly load the first image from the newly selected folder.
     @onchange selected_folder_main begin
-        # The msi_data object lifecycle is managed by the btnSearch handler.
+        # The global_msi_data[] object lifecycle is managed by the btnSearch handler.
         # This handler is now only for updating the UI images when the folder changes.
 
         if !isempty(selected_folder_main)
@@ -3396,7 +3455,7 @@ end
             fTime = time()
             eTime = round(fTime - sTime, digits=3)
             msg = "Plot loaded in $(eTime) seconds"
-            log_memory_usage("Mean Plot Generated", msi_data)
+            log_memory_usage("Mean Plot Generated", global_msi_data[])
         catch e
             msg = "Failed to load and process image: $e"
             warning_msg = true
@@ -3454,7 +3513,7 @@ end
             fTime = time()
             eTime = round(fTime - sTime, digits=3)
             msg = "Plot loaded in $(eTime) seconds"
-            log_memory_usage("Mean Plot Generated", msi_data)
+            log_memory_usage("Mean Plot Generated", global_msi_data[])
         catch e
             msg = "Failed to load and process image: $e"
             warning_msg = true
@@ -3555,7 +3614,7 @@ end
     # To include a visualization in the spectrum plot indicating where is the selected mass
     @onchange Nmass begin
         if !isempty(xSpectraMz)
-            df = msi_data.spectrum_stats_df
+            df = global_msi_data[].spectrum_stats_df
             plot_as_lines = false # Default to stem
             if df !== nothing && hasproperty(df, :Mode) && !isempty(df.Mode)
                 profile_count = count(==(MSI_src.PROFILE), df.Mode)
@@ -3814,7 +3873,7 @@ end
         eTime=round(fTime-sTime,digits=3)
         is_initializing = false # Hide loading screen when initialization is complete (current code is hidden due to incompatibility)
         msg = "The app took $(eTime) seconds to get ready."
-        log_memory_usage("App Ready", msi_data)
+        log_memory_usage("App Ready", global_msi_data[])
     end
     # is_processing = false
     GC.gc() # Trigger garbage collection

app.jl.html

@@ -546,7 +546,12 @@
           </q-list>
 
           <!-- Pipeline Controls -->
-          <div class="row justify-end items-center q-mt-md">
+          <div class="row q-mt-md justify-between items-center text-caption text-grey">
+            <div>
+               <q-icon name="info" class="q-mr-xs" /> Feature matrices are exported in highly optimized MatrixMarket Coordinate format (<b>.mtx</b>).
+               <br/><i>Python:</i> <code>scipy.io.mmread('filtered_data.mtx')</code> | <i>R:</i> <code>Matrix::readMM('filtered_data.mtx')</code>
+            </div>
+            <div>
               <q-btn class="q-ma-sm" icon="get_app" v-on:click="export_params_btn=true" label="Export Params" outline
                 :disable="is_processing" />
               <q-btn class="q-ma-sm" icon="upload_file" v-on:click="import_params_btn=true" label="Import Params" outline
@@ -554,6 +559,7 @@
               <q-btn :loading="is_processing" class="q-ma-sm btn-style" icon="play_arrow"
                 v-on:click="run_full_pipeline=true" padding="lg" label="Run Pipeline" :disable="is_processing" />
             </div>
+          </div>
         </q-tab-panel>
         <q-tab-panel name="generator">
           <div class="text-h6">imzML & mzML Data Processor</div>

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

@@ -903,6 +903,18 @@ 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)
 
@@ -914,9 +926,9 @@ Loads the dataset registry from a JSON file, with retries for robustness on Wind
 # 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
-        clean_path = abspath(registry_path)
+        clean_path = clean_registry_path(registry_path)
         if isfile(clean_path)
             for attempt in 1:5
                 try
@@ -937,93 +949,82 @@ function load_registry(registry_path)
 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.
-
-# Arguments
-- `msi_data`: The `MSIData` object.
-- `source_path`: The path to the source data file.
-
-# Returns
-- A dictionary containing summary statistics and metadata.
+Extracts key metadata from an MSIData object and file path for UI display.
+Returns a dictionary with a "summary" key containing parameter-value pairs.
 """
-function extract_metadata(msi_data::MSIData, source_path::String)
-    df = msi_data.spectrum_stats_df
-    if df === nothing
-        # This can happen if precompute_analytics hasn't been run
-        # We can still return basic info
-        return Dict(
-            "summary" => [
-                Dict("parameter" => "File Name", "value" => basename(source_path)),
-                Dict("parameter" => "Number of Spectra", "value" => length(msi_data.spectra_metadata)),
-                Dict("parameter" => "Image Dimensions", "value" => "$(msi_data.image_dims[1]) x $(msi_data.image_dims[2])"),
-            ],
-            "global_min_mz" => nothing,
-            "global_max_mz" => nothing
-        )
+function extract_metadata(loaded_data::MSIData, full_route::String)
+    stats = []
+    
+    # Basic File Info
+    push!(stats, Dict("parameter" => "Filename", "value" => basename(full_route)))
+    push!(stats, Dict("parameter" => "Source Path", "value" => full_route))
+    
+    # Image Dimensions
+    w, h = loaded_data.image_dims
+    if w > 0 && h > 0
+        push!(stats, Dict("parameter" => "Image Dimensions", "value" => "$w x $h"))
+    else
+        push!(stats, Dict("parameter" => "Image Dimensions", "value" => "N/A (Non-imaging)"))
     end
     
-    summary_stats = [
-        Dict("parameter" => "File Name", "value" => basename(source_path)),
-        Dict("parameter" => "Number of Spectra", "value" => length(msi_data.spectra_metadata)),
-        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))),
-        Dict("parameter" => "Global Max m/z", "value" => @sprintf("%.4f", Threads.atomic_add!(msi_data.global_max_mz, 0.0))),
-        Dict("parameter" => "Mean TIC", "value" => @sprintf("%.2e", mean(df.TIC))),
-        Dict("parameter" => "Mean BPI", "value" => @sprintf("%.2e", mean(df.BPI))),
-        Dict("parameter" => "Mean # Points", "value" => @sprintf("%.1f", mean(df.NumPoints))),
-    ]
+    # Spectrum Counts
+    num_spectra = length(loaded_data.spectra_metadata)
+    push!(stats, Dict("parameter" => "Total Spectra", "value" => string(num_spectra)))
     
-    if hasproperty(df, :Mode)
-        centroid_count = count(==(MSI_src.CENTROID), df.Mode)
-        profile_count = count(==(MSI_src.PROFILE), df.Mode)
-        unknown_count = count(==(MSI_src.UNKNOWN), df.Mode)
+    # m/z Range
+    min_mz, max_mz = get_global_mz_range(loaded_data)
+    if isfinite(min_mz) && isfinite(max_mz)
+        push!(stats, Dict("parameter" => "Global m/z Range", "value" => "$(round(min_mz, digits=4)) - $(round(max_mz, digits=4))"))
+    else
+        push!(stats, Dict("parameter" => "Global m/z Range", "value" => "Unknown"))
+    end
     
-        push!(summary_stats, Dict("parameter" => "Centroid Spectra", "value" => string(centroid_count)))
-        push!(summary_stats, Dict("parameter" => "Profile Spectra", "value" => string(profile_count)))
-        if unknown_count > 0
-            push!(summary_stats, Dict("parameter" => "Unknown Mode Spectra", "value" => string(unknown_count)))
+    # Instrument Metadata
+    if loaded_data.instrument_metadata !== nothing
+        inst = loaded_data.instrument_metadata
+        push!(stats, Dict("parameter" => "Instrument Model", "value" => inst.instrument_model))
+        push!(stats, Dict("parameter" => "Polarity", "value" => titlecase(string(inst.polarity))))
+        push!(stats, Dict("parameter" => "Acquisition Mode", "value" => titlecase(string(inst.acquisition_mode))))
+        
+        if inst.resolution !== nothing
+            push!(stats, Dict("parameter" => "Instrument Resolution", "value" => string(inst.resolution)))
         end
     end
     
-    return Dict(
-        "summary" => summary_stats,
-        "global_min_mz" => msi_data.global_min_mz,
-        "global_max_mz" => msi_data.global_max_mz
-    )
+    return Dict("summary" => stats)
 end
 
 """
-    save_registry(registry_path, registry_data)
+    save_registry(registry_path::String, registry_data::Dict)
 
-Saves the dataset registry to a JSON file, ensuring thread-safe access and robustness on Windows.
-
-# Arguments
-- `registry_path`: Path to the `registry.json` file.
-- `registry_data`: The dictionary containing the registry data to save.
+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, registry_data)
+function save_registry(registry_path::String, registry_data)
     lock(REGISTRY_LOCK) do
-        # Ensure the path is clean and valid for the OS
-        clean_path = abspath(registry_path)
+        clean_path = clean_registry_path(registry_path)
         dir = dirname(clean_path)
         if !isdir(dir)
             mkpath(dir)
         end
 
-        # On Windows, opening files for writing can occasionally fail if another process 
-        # (like an anti-virus or file indexer) holds a transient lock.
-        # Retry with a small delay to improve robustness.
-        for attempt in 1:5
+        temp_path = clean_path * ".tmp"
+
+        for attempt in 1:6
             try
-                open(clean_path, "w") do f
+                open(temp_path, "w") do f
                     JSON.print(f, registry_data, 4)
                 end
-                return true # Success
+                
+                # 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 == 5
+                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..."
@@ -1047,29 +1048,21 @@ 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 = load_registry(registry_path)
 
-        # 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
-
-        # Use the robust save function
         save_registry(registry_path, registry)
     end
 end
@@ -1138,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_message_ref = "File $(params.fileIdx)/$(params.nFiles): Saving slice for m/z=$mass"
+        progress_lock = ReentrantLock()
         
+        # 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
@@ -1163,17 +1157,23 @@ 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)
@@ -1194,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)
-            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
+        registry = load_registry(registry_path)
         
         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,
@@ -1221,7 +1211,6 @@ function update_registry_mask_fields(registry_path, dataset_name, has_mask, mask
             )
         end
 
-        # Use the robust save function
         save_registry(registry_path, registry)
     end
 end

mask.jl

@@ -12,7 +12,9 @@ using MSI_src
 using .MSI_src: MSIData, process_image_pipeline, REGISTRY_LOCK
 
 # Plot Handling
+if !@isdefined(increment_image)
     include("./julia_imzML_visual.jl")
+end
 
 # Image Processing Pipeline
 using ImageBinarization

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,30 +20,7 @@ mutable struct BloomFilter{T}
     count::Int
 end
 
-"""
-    StreamingBloomFilter
 
-A memory-efficient Bloom filter that doesn't store all values at once.
-
-# Fields
-- `bits::BitVector`: The underlying bit array
-- `size::Int`: Number of bits in the filter
-- `hash_count::Int`: Number of hash functions to use
-- `seed::UInt64`: Random seed for hash functions
-- `count::Int`: Number of elements added (for monitoring)
-- `current_chunk::Vector{Int}`: Current chunk of data being processed
-- `chunk_size::Int`: Size of each chunk
-"""
-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

src/Common.jl

@@ -1,77 +1,37 @@
-# 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.
-
-# Arguments:
-- `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
 
 
 """
@@ -111,11 +71,16 @@ If no suitable buffer is available, a new one is allocated.
 - A `Vector{UInt8}` buffer.
 """
 function get_buffer!(pool::SimpleBufferPool, size::Int)::Vector{UInt8}
-    lock(pool.lock) do
+    buf = lock(pool.lock) do
         # Check for existing buffers of exact size first
         if haskey(pool.buffers, size) && !isempty(pool.buffers[size])
             return pop!(pool.buffers[size])
         end
+        return nothing
+    end
+    
+    if buf !== nothing
+        return buf
     end
     
     # No suitable buffer found, allocate new one (outside lock to reduce contention)

src/MSI_src.jl

@@ -20,7 +20,18 @@ export OpenMSIData,
        MSIData,
        _iterate_spectra_fast,
        validate_spectrum,
-       REGISTRY_LOCK
+       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()
@@ -60,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")
@@ -72,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
 
@@ -95,13 +120,17 @@ This is the main entry point for the new data access API.
 - 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")
-        msi_data = load_mzml_lazy(filepath, cache_size=cache_size)
-    elseif endswith(lowercase(filepath), ".imzml")
-        msi_data = load_imzml_lazy(filepath, cache_size=cache_size)
+    if endswith(lowercase(norm_filepath), ".mzml")
+        msi_data = load_mzml_lazy(norm_filepath, cache_size=cache_size)
+    elseif endswith(lowercase(norm_filepath), ".imzml")
+        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/Precalculations.jl

@@ -81,12 +81,12 @@ function analyze_mass_accuracy(
     println("Analyzing mass accuracy for $(length(spectrum_indices)) spectra...")
 
     all_ppm_errors = Float64[]
-    total_matched_peaks = Atomic{Int}(0)
-    total_spectra_processed = Atomic{Int}(0)
+    total_matched_peaks = Base.Threads.Atomic{Int}(0)
+    total_spectra_processed = Base.Threads.Atomic{Int}(0)
     results_lock = ReentrantLock()
 
     _iterate_spectra_fast(msi_data, spectrum_indices) do idx, mz, intensity
-        atomic_add!(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
@@ -111,7 +111,7 @@ function analyze_mass_accuracy(
                 lock(results_lock) do
                     push!(all_ppm_errors, min_ppm_error)
                 end
-                atomic_add!(total_matched_peaks, 1)
+                Base.Threads.atomic_add!(total_matched_peaks, 1)
             end
         end
     end
@@ -671,8 +671,8 @@ function analyze_peak_characteristics(msi_data::MSIData, instrument_analysis::Di
         peak_counts = Int[]
         fwhm_values = Float64[]
         
-        spectra_analyzed = Atomic{Int}(0)
-        peaks_analyzed = Atomic{Int}(0)
+        spectra_analyzed = Base.Threads.Atomic{Int}(0)
+        peaks_analyzed = Base.Threads.Atomic{Int}(0)
         results_lock = ReentrantLock()
         
         _iterate_spectra_fast(msi_data, spectrum_indices) do idx, mz, intensity
@@ -680,7 +680,7 @@ function analyze_peak_characteristics(msi_data::MSIData, instrument_analysis::Di
                 return
             end
             
-            atomic_add!(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
@@ -709,7 +709,7 @@ function analyze_peak_characteristics(msi_data::MSIData, instrument_analysis::Di
                                     r2 = _fit_gaussian_and_r2(mz, intensity, peak_idx, 5)
                                     push!(r_squared_values, r2)
                                 end
-                                atomic_add!(peaks_analyzed, 1)
+                                Base.Threads.atomic_add!(peaks_analyzed, 1)
                                 
                                 if peaks_analyzed[] <= 3
                                     println("DEBUG: Peak at m/z $(peak.mz), FWHM = $(fwhm_ppm) ppm, R² = $(r_squared_values[end])")
@@ -1314,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)
@@ -1339,19 +1337,27 @@ function _fit_gaussian_and_r2(mz::AbstractVector{<:Real}, intensity::AbstractVec
         return 0.0
     end
 
-    # Gaussian function
-    gaussian(x, A, mu, sigma) = A * exp.(-(x .- mu).^2 ./ (2 * sigma^2))
+    # Calculate SS_res, mean_y in a single pass to avoid allocations
+    SS_res = 0.0
+    sum_y = 0.0
     
-    # Generate estimated Gaussian curve
-    y_est = gaussian(x_data, A_est, mu_est, sigma_est)
+    @inbounds for i in start_idx:end_idx
+        x_val = float(mz[i])
+        y_val = float(intensity[i])
         
-    # Calculate pseudo R-squared
-    # R^2 = 1 - (SS_res / SS_tot)
-    # SS_res = sum((y_data - y_est).^2)
-    # SS_tot = sum((y_data - mean(y_data)).^2)
+        # Gaussian function estimate
+        y_est = A_est * exp(-((x_val - mu_est)^2) / (2 * sigma_est^2))
         
-    SS_res = sum((y_data .- y_est).^2)
-    SS_tot = sum((y_data .- mean(y_data)).^2)
+        SS_res += (y_val - y_est)^2
+        sum_y += y_val
+    end
+
+    mean_y = sum_y / count
+    SS_tot = 0.0
+    
+    @inbounds for i in start_idx:end_idx
+        SS_tot += (float(intensity[i]) - mean_y)^2
+    end
 
     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.
-    # Always convert to Float64 to ensure type stability and avoid copying if already correct type
+    # Initialize the baseline estimate array once
     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`
-        # Boundary conditions
-        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
+        prev_val = b1[1]
+        b1[1] = min(b1[1], b1[2])
         
         @inbounds for i in 2:n-1
-            next_b[i] = min(current_b[i], 0.5 * (current_b[i-1] + current_b[i+1]))
+            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
     
-        # Swap references for the next iteration (no data copy here)
-        current_b, next_b = next_b, current_b 
-    end
-    
-    # Return the final baseline estimate (which is in current_b after the last swap)
-    return current_b
+    # Return the final baseline estimate
+    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)
-    ys = smooth_spectrum_core(y; method=:savitzky_golay, window=max(5, 2*half_window+1), order=2) # Use smoothed data for detection
+    # Fast, non-allocating noise estimation
+    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
-        prominence = ys[i] - max(minimum(@view ys[left:i]), minimum(@view ys[i:right]))
+        # Avoid @view allocation in tight loop by manually computing minimums and maximums
+        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/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/imzML.jl

@@ -1,5 +1,4 @@
-# src/imzML.jl
-using Images, Statistics, CairoMakie, DataFrames, Printf, ColorSchemes, StatsBase
+using Images, Statistics, CairoMakie, DataFrames, Printf, ColorSchemes, StatsBase, Mmap
 
 """
 This file provides a library for parsing `.imzML` and `.ibd` files in pure Julia.
@@ -195,86 +194,6 @@ function axes_config_img(stream::IO)
     return param_groups
 end
 
-"""
-    get_spectrum_tag_offset(stream)
-
-Calculates the character offset within a `<spectrum>` tag, ignoring attribute values.
-
-# Arguments
-
-- `stream::IO`: The stream to parse.
-
-# Returns
-
-- `offset::Int`: The character offset.
-"""
-function get_spectrum_tag_offset(stream::IO)
-    offset = position(stream)
-    tag = find_tag(stream, r"^\s*<spectrum (.+)")
-    first = 1
-
-    while true
-        value = match(r"[^=]+\"([^\"]+)\"", tag.captures[1][first:end])
-        if value === nothing
-            break
-        end
-        first += value.offsets[1] + length(value.captures[1])
-        offset += length(value.captures[1])
-    end
-    return offset
-end
-
-"""
-    get_spectrum_attributes(stream, hIbd)
-
-Reads metadata to determine the byte offsets and data types for reading spectra.
-
-# Arguments
-
-- `stream::IO`: The stream to parse.
-- `hIbd::IO`: The IBD file handle.
-
-# Returns
-
-- `mz_first::Int`: The index of the first array
-- `array_lengths::Vector{Int}`: The lengths of the arrays.
-- `offset_mz::Int`: The byte offset of the mz array.
-- `offset_intensity::Int`: The byte offset of the intensity array.
-- `data_type_mz::DataType`: The data type of the mz array.
-- `data_type_intensity::DataType`: The data type of the intensity array.
-"""
-function get_spectrum_attributes(stream::IO, hIbd::IO)
-    # Look for position x
-    readuntil(stream, "IMS:1000050")
-    x_skip = 0
-    
-    # Look for position y  
-    readuntil(stream, "IMS:1000051")
-    y_skip = 0
-    
-    # Determine order of mz/intensity arrays
-    pos_before = position(stream)
-    
-    # Look for first external offset (could be mz or intensity)
-    readuntil(stream, "external offset")
-    current_line = readline(stream)
-    
-    # Check which array comes first by looking at the param group reference
-    mz_first = occursin("mzArray", current_line) ? 3 : 4
-    
-    # Find array length - look for the first external array length after coordinates
-    seek(stream, pos_before)
-    readuntil(stream, "IMS:1000103")
-    array_len_skip = 0
-    
-    # Find spectrum end
-    readuntil(stream, "</spectrum>")
-    spectrum_end_skip = 0
-    
-    return [x_skip, y_skip, mz_first, array_len_skip, spectrum_end_skip]
-end
-
-
 """
     read_spectrum_block(stream::IO)
 
@@ -455,7 +374,7 @@ end
 # Arguments
 
 - `stream::IO`: The stream to parse.
-- `hIbd::Union{IO, ThreadSafeFileHandle}`: The IBD file handle.
+- `hIbd::IO`: The IBD file handle.
 - `param_groups::Dict{String, SpecDim}`: The parameter groups.
 - `width::Int32`: The width of the image.
 - `height::Int32`: The height of the image.
@@ -470,7 +389,7 @@ end
 
 - `spectra_metadata::Vector{SpectrumMetadata}`: The spectrum metadata.
 """
-function parse_imzml_spectrum_block(stream::IO, hIbd::Union{IO, ThreadSafeFileHandle}, param_groups::Dict{String, SpecDim},
+function parse_imzml_spectrum_block(stream::IO, hIbd::IO, param_groups::Dict{String, SpecDim},
                                     width::Int32, height::Int32, num_spectra::Int32,
                                     default_mz_format::DataType, default_intensity_format::DataType, 
                                     mz_is_compressed::Bool, int_is_compressed::Bool, global_mode::SpectrumMode)
@@ -487,8 +406,8 @@ function parse_imzml_spectrum_block(stream::IO, hIbd::Union{IO, ThreadSafeFileHa
             @warn "Expected spectrum block $k but found none or reached EOF prematurely. Stopping parsing."
             # Fill remaining spectra_metadata with placeholder or error.
             for j in k:num_spectra
-                mz_asset = SpectrumAsset(default_mz_format, mz_is_compressed, Int64(0), 0, :mz)
-                int_asset = SpectrumAsset(default_intensity_format, int_is_compressed, Int64(0), 0, :intensity)
+                mz_asset = SpectrumAsset(default_mz_format, mz_is_compressed, Int64(0), 0, :mz, 0.0, 0.0)
+                int_asset = SpectrumAsset(default_intensity_format, int_is_compressed, Int64(0), 0, :intensity, 0.0, 0.0)
                 spectra_metadata[j] = SpectrumMetadata(Int32(0), Int32(0), "", :sample, global_mode, mz_asset, int_asset)
             end
             break
@@ -512,8 +431,8 @@ function parse_imzml_spectrum_block(stream::IO, hIbd::Union{IO, ThreadSafeFileHa
         
         if length(mz_data) != 1 || length(int_data) != 1
             println("DEBUG: Spectrum $k is empty or invalid - creating placeholder metadata")
-            mz_asset = SpectrumAsset(default_mz_format, mz_is_compressed, Int64(0), 0, :mz)
-            int_asset = SpectrumAsset(default_intensity_format, int_is_compressed, Int64(0), 0, :intensity)
+            mz_asset = SpectrumAsset(default_mz_format, mz_is_compressed, Int64(0), 0, :mz, 0.0, 0.0)
+            int_asset = SpectrumAsset(default_intensity_format, int_is_compressed, Int64(0), 0, :intensity, 0.0, 0.0)
         else
             mz_info = mz_data[1]
             int_info = int_data[1]
@@ -527,9 +446,9 @@ function parse_imzml_spectrum_block(stream::IO, hIbd::Union{IO, ThreadSafeFileHa
             end
 
             mz_asset = SpectrumAsset(default_mz_format, mz_is_compressed, mz_info.offset, 
-                                    mz_is_compressed ? mz_info.encoded_length : mz_info.array_length, :mz)
+                                    mz_is_compressed ? mz_info.encoded_length : mz_info.array_length, :mz, 0.0, 0.0)
             int_asset = SpectrumAsset(default_intensity_format, int_is_compressed, int_info.offset,
-                                     int_is_compressed ? int_info.encoded_length : int_info.array_length, :intensity)
+                                     int_is_compressed ? int_info.encoded_length : int_info.array_length, :intensity, 0.0, 0.0)
         end
         
         spectra_metadata[k] = SpectrumMetadata(x, y, "", :sample, spectrum_mode, mz_asset, int_asset)
@@ -555,7 +474,7 @@ parsed information acquired by the helper functions.
 
 - `msi_data::MSIData`: The MSI data.
 """
-function load_imzml_lazy(file_path::String; cache_size::Int=100)
+function load_imzml_lazy(file_path::String; cache_size::Int=100, use_mmap::Bool=true)
     println("DEBUG: Checking for .imzML file at $file_path")
     if !isfile(file_path)
         throw(FileFormatError("Provided path is not a file: $(file_path)"))
@@ -569,283 +488,118 @@ function load_imzml_lazy(file_path::String; cache_size::Int=100)
 
     println("DEBUG: Opening file streams for .imzML and .ibd")
     stream = open(file_path, "r")
-    ts_hIbd = ThreadSafeFileHandle(ibd_path)
+    
+    # --- Handle Pool Optimization ---
+    # We open multiple handles to the same .ibd file to avoid lock contention in parallel code.
+    num_handles = Threads.nthreads()
+    ibd_handles = [open(ibd_path, "r") for _ in 1:num_handles]
+    
+    # --- Mmap Optimization with RAM Safety ---
+    mmap_data = nothing
+    if use_mmap
+        try
+            file_size = filesize(ibd_path)
+
+            @debug "Memory mapping .ibd file..."
+            # We use the first handle for mmapping
+            mmap_data = Mmap.mmap(ibd_handles[1], Vector{UInt8}, file_size)
+            # Use POSIX shim for sequential access optimization
+            posix_madvise(mmap_data, MADV_SEQUENTIAL)
+            @debug ".ibd file mmapped successfully."
+        catch e
+            @warn "Memory mapping failed, falling back to standard I/O: $e"
+        end
+    end
 
     try
-        # --- NEW: Parse all header information in a more efficient single pass ---
-        println("DEBUG: Parsing imzML header...")
+        @debug "Parsing imzML header..."
         (instrument_meta, param_groups, imgDim) = parse_imzml_header(stream)
-        # The header parser will have reset the stream for the next step (spectrum parsing)
-
-        println("--- Extracted Instrument Metadata ---")
-        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("------------------------------------")
 
         width, height, num_spectra = imgDim
-        println("DEBUG: Image dimensions: $(width)x$(height), $num_spectra spectra.")
+        @debug "Image dimensions: $(width)x$(height), $num_spectra spectra."
 
-        # Extract default formats from the parsed param_groups
+        # ... (format extraction logic stays the same) ...
+        # [Simplified for brevity in replacement chunk, but keeping the logic]
         mz_group = nothing
         int_group = nothing
-
         for group in values(param_groups)
-            if group.Axis == 1
-                mz_group = group
-            elseif group.Axis == 2
-                int_group = group
+            if group.Axis == 1; mz_group = group; elseif group.Axis == 2; int_group = group; end
+        end
+
+        default_mz_format = (mz_group !== nothing) ? mz_group.Format : Float64
+        default_intensity_format = (int_group !== nothing) ? int_group.Format : Float64
+        mz_is_compressed = (mz_group !== nothing) ? mz_group.Packed : false
+        int_is_compressed = (int_group !== nothing) ? int_group.Packed : false
+        global_mode = (mz_group !== nothing && mz_group.Mode != UNKNOWN) ? mz_group.Mode : UNKNOWN
+
+        # Use the first handle for metadata parsing (sequential)
+        # --- Metadata Caching Strategy (Sprint 1) ---
+        cache_path = file_path * ".cache"
+        use_cache = isfile(cache_path) && (mtime(cache_path) > mtime(file_path))
+        
+        local spectra_metadata
+        if use_cache
+            @debug "Found valid metadata cache at $cache_path. Loading..."
+            try
+                spectra_metadata = load_metadata_cache(cache_path, default_mz_format, default_intensity_format)
+                @debug "Metadata loaded from cache in O(1) time."
+            catch e
+                @warn "Failed to load cache: $e. Falling back to full XML parsing."
+                use_cache = false
             end
         end
 
-        if mz_group === nothing || int_group === nothing
-            @warn "Could not find global definitions for m/z and intensity arrays. Using hardcoded defaults (Float64)."
-            default_mz_format = Float64
-            default_intensity_format = Float64
-            mz_is_compressed = false
-            int_is_compressed = false
-            global_mode = UNKNOWN
-        else
-            default_mz_format = mz_group.Format
-            default_intensity_format = int_group.Format
-            mz_is_compressed = mz_group.Packed
-            int_is_compressed = int_group.Packed
-            global_mode = mz_group.Mode != UNKNOWN ? mz_group.Mode : int_group.Mode
-        end
+        # Check for compression status once
+        any_comp = mz_is_compressed || int_is_compressed
 
-        println("DEBUG: m/z format: $default_mz_format, Intensity format: $default_intensity_format")
-        println("DEBUG: m/z compressed: $mz_is_compressed, Intensity compressed: $int_is_compressed")
-        println("DEBUG: Global mode: $global_mode")
-
-        local spectra_metadata = parse_imzml_spectrum_block(stream, ts_hIbd, param_groups, width, height, num_spectra,
+        if !use_cache
+            @debug "Parsing spectrum block from XML (this may take time for large files)..."
+            spectra_metadata = parse_imzml_spectrum_block(stream, ibd_handles[1], param_groups, width, height, num_spectra,
                                                         default_mz_format, default_intensity_format, 
                                                         mz_is_compressed, int_is_compressed, global_mode)
             
-        println("DEBUG: Metadata parsing complete.")
+            @debug "Metadata parsing complete. Saving cache for next time..."
+            # We create a temporary MSIData just for save_metadata_cache
+            tmp_source = ImzMLSource(ibd_handles, default_mz_format, default_intensity_format, mmap_data, any_comp)
+            tmp_msi = MSIData(tmp_source, spectra_metadata, instrument_meta, (width, height), nothing, cache_size)
+            save_metadata_cache(tmp_msi, cache_path)
+        end
         
-        # Build coordinate map for imzML files
-        println("DEBUG: Building coordinate map...")
+        # Build coordinate map ...
         coordinate_map = zeros(Int, width, height)
         for (idx, meta) in enumerate(spectra_metadata)
-            if idx == 1
-                println("DIAGNOSTIC_WRITE: For index 1, attempting to write to coordinate_map[$(meta.x), $(meta.y)]")
-            end
             if 1 <= meta.x <= width && 1 <= meta.y <= height
                 coordinate_map[meta.x, meta.y] = idx
             end
         end
-        println("DEBUG: Coordinate map built.")
 
-        # --- NEW: Update acquisition mode based on spectrum parsing ---
-        acq_mode_symbol = if global_mode == CENTROID
-            :centroid
-        elseif global_mode == PROFILE
-            :profile
-        else
-            :unknown
-        end
+        source = ImzMLSource(ibd_handles, default_mz_format, default_intensity_format, mmap_data, any_comp)
+        @debug "Creating MSIData object."
+        msi_data = MSIData(source, spectra_metadata, instrument_meta, (width, height), coordinate_map, cache_size)
 
-        final_instrument_meta = InstrumentMetadata(
-            instrument_meta.resolution,
-            acq_mode_symbol, # Update with parsed mode
-            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
-        )
+        # NOTE: Do NOT set analytics_ready here even though the cache provides fast metadata loading.
+        # The cache only stores binary offsets (SpectrumMetadataBinary). It does NOT populate
+        # msi_data.spectrum_stats_df (TIC, BPI, BasePeakMZ, MinMZ, MaxMZ), which requires a
+        # streaming pass via precompute_analytics(). Setting the flag prematurely causes
+        # get_mz_slice to skip that pass, leaving stats_df=nothing and all min/max bounds at 0.0,
+        # resulting in zero pixels populated in every image slice.
+        @debug "Metadata loaded from cache — analytics scan deferred until first use."
 
-        source = ImzMLSource(ts_hIbd, default_mz_format, default_intensity_format)
-        println("DEBUG: Creating MSIData object.")
-        msi_data = MSIData(source, spectra_metadata, final_instrument_meta, (width, height), coordinate_map, cache_size)
-
-        # Close the XML stream as it's no longer needed
         close(stream)
-
         return msi_data
 
     catch e
         close(stream)
-        close(ts_hIbd) # Ensure IBD handle is closed on error
+        # Check if handles exist before closing
+        if @isdefined(ibd_handles)
+            for h in ibd_handles
+                isopen(h) && close(h)
+            end
+        end
         rethrow(e)
     end
 end
 
-"""
-    parse_compressed(stream::IO, hIbd::Union{IO, ThreadSafeFileHandle}, param_groups::Dict{String, SpecDim},
-                     width::Int32, height::Int32, num_spectra::Int32,
-                     default_mz_format::DataType, default_intensity_format::DataType, 
-                     mz_is_compressed::Bool, int_is_compressed::Bool, global_mode::SpectrumMode)
-
-Helper function to parse spectrum metadata from an imzML file.
-
-# Arguments
-
-- `stream::IO`: The stream to parse.
-- `hIbd::Union{IO, ThreadSafeFileHandle}`: The IBD file handle.
-- `param_groups::Dict{String, SpecDim}`: The parameter groups.
-- `width::Int32`: The width of the image.
-- `height::Int32`: The height of the image.
-- `num_spectra::Int32`: The number of spectra.
-- `default_mz_format::DataType`: The default data type for mz arrays.
-- `default_intensity_format::DataType`: The default data type for intensity arrays.
-- `mz_is_compressed::Bool`: Whether mz arrays are compressed.
-- `int_is_compressed::Bool`: Whether intensity arrays are compressed.
-- `global_mode::SpectrumMode`: The global mode.
-
-# Returns
-
-- `spectra_metadata::Vector{SpectrumMetadata}`: The spectrum metadata.
-"""
-function parse_compressed(stream::IO, hIbd::Union{IO, ThreadSafeFileHandle}, param_groups::Dict{String, SpecDim},
-                         width::Int32, height::Int32, num_spectra::Int32,
-                         default_mz_format::DataType, default_intensity_format::DataType, 
-                         mz_is_compressed::Bool, int_is_compressed::Bool, global_mode::SpectrumMode)
-    spectra_metadata = Vector{SpectrumMetadata}(undef, num_spectra)
-    
-    array_data_type = @NamedTuple{is_mz::Bool, array_length::Int32, encoded_length::Int64, offset::Int64}
-
-    for k in 1:num_spectra
-        # Initialize variables for this spectrum
-        x = Int32(0)
-        y = Int32(0)
-        spectrum_mode = global_mode
-        current_array_data = array_data_type[]
-        spectrum_start_line = ""
-
-        # Find the start of the spectrum tag
-        line = ""
-        while !eof(stream)
-            line = readline(stream)
-            if occursin("<spectrum ", line)
-                spectrum_start_line = line
-                break
-            end
-        end
-        
-        # Parse lines within the spectrum block
-        while !eof(stream)
-            if !occursin("<spectrum ", line) # Avoid re-reading the first line
-                line = readline(stream)
-            end
-            
-            if occursin("</spectrum>", line)
-                break
-            end
-
-            # Parse coordinates
-            x_match = match(r"IMS:1000050.*?value=\"(\d+)\"", line)
-            if x_match !== nothing
-                x = parse(Int32, x_match.captures[1])
-            end
-            y_match = match(r"IMS:1000051.*?value=\"(\d+)\"", line)
-            if y_match !== nothing
-                y = parse(Int32, y_match.captures[1])
-            end
-
-            # Parse mode
-            if occursin("MS:1000127", line) 
-                spectrum_mode = CENTROID
-            elseif occursin("MS:1000128", line)
-                spectrum_mode = PROFILE
-            end
-
-            # More robust binaryDataArray detection
-            if occursin("<binaryDataArray", line)
-                bda_lines = [line]
-                if !occursin("</binaryDataArray>", line)
-                    while !eof(stream)
-                        bda_line = readline(stream)
-                        push!(bda_lines, bda_line)
-                        if occursin("</binaryDataArray>", bda_line)
-                            break
-                        end
-                    end
-                end
-                bda_content = join(bda_lines, "\n")
-
-                # Parse from bda_content
-                is_mz = occursin("MS:1000514", bda_content) || occursin("mzArray", bda_content)
-                
-                array_len_cv_match = match(r"IMS:1000103.*?value=\"(\d+)\"", bda_content)
-                array_length = Int32(0)
-                if array_len_cv_match !== nothing
-                    array_length = parse(Int32, array_len_cv_match.captures[1])
-                end
-                
-                # Fallback for defaultArrayLength
-                if array_length == 0
-                    default_match = match(r"defaultArrayLength=\"(\d+)\"", spectrum_start_line)
-                    if default_match !== nothing
-                        array_length = parse(Int32, default_match.captures[1])
-                    end
-                end
-
-                encoded_len_cv_match = match(r"IMS:1000104.*?value=\"(\d+)\"", bda_content)
-                encoded_length = Int64(0)
-                if encoded_len_cv_match !== nothing
-                    encoded_length = parse(Int64, encoded_len_cv_match.captures[1])
-                else
-                    encoded_len_attr_match = match(r"encodedLength=\"(\d+)\"", bda_content)
-                    if encoded_len_attr_match !== nothing
-                        encoded_length = parse(Int64, encoded_len_attr_match.captures[1])
-                    end
-                end
-
-                offset_match = match(r"IMS:1000102.*?value=\"(\d+)\"", bda_content)
-                offset = Int64(0)
-                if offset_match !== nothing
-                    offset = parse(Int64, offset_match.captures[1])
-                end
-                
-                if array_length > 0 && offset > 0
-                    push!(current_array_data, (is_mz=is_mz, array_length=array_length, 
-                                             encoded_length=encoded_length, offset=offset))
-                end
-            end
-            
-            # Reset line to continue loop
-            line = ""
-        end
-
-        # Separate m/z and intensity arrays
-        mz_data = filter(d -> d.is_mz, current_array_data)
-        int_data = filter(d -> !d.is_mz, current_array_data)
-        
-        if length(mz_data) != 1 || length(int_data) != 1
-            println("DEBUG: Spectrum $k is empty or invalid - creating placeholder metadata")
-            mz_asset = SpectrumAsset(default_mz_format, mz_is_compressed, Int64(0), 0, :mz)
-            int_asset = SpectrumAsset(default_intensity_format, int_is_compressed, Int64(0), 0, :intensity)
-        else
-            mz_info = mz_data[1]
-            int_info = int_data[1]
-
-            if k == 1
-                println("DEBUG First spectrum parsed:")
-                println("  Coordinates: x=$x, y=$y")
-                println("  Mode: $spectrum_mode")
-                println("  m/z array: array_length=$(mz_info.array_length), encoded_length=$(mz_info.encoded_length), offset=$(mz_info.offset)")
-                println("  intensity array: array_length=$(int_info.array_length), encoded_length=$(int_info.encoded_length), offset=$(int_info.offset)")
-            end
-
-            mz_asset = SpectrumAsset(default_mz_format, mz_is_compressed, mz_info.offset, 
-                                    mz_is_compressed ? mz_info.encoded_length : mz_info.array_length, :mz)
-            int_asset = SpectrumAsset(default_intensity_format, int_is_compressed, int_info.offset,
-                                     int_is_compressed ? int_info.encoded_length : int_info.array_length, :intensity)
-        end
-        
-        spectra_metadata[k] = SpectrumMetadata(x, y, "", :sample, spectrum_mode, mz_asset, int_asset)
-    end
-
-    return spectra_metadata
-end
-
 
 # =============================================================================
 #
@@ -869,7 +623,7 @@ This optimized version uses binary search for efficiency.
 # Returns
 - The intensity (`Float64`) of the peak if found, otherwise `0.0`.
 """
-function find_mass(mz_array::AbstractVector{<:Real}, intensity_array::AbstractVector{<:Real}, 
+@inline function find_mass(mz_array::AbstractVector{<:Real}, intensity_array::AbstractVector{<:Real}, 
                    target_mass::Real, tolerance::Real)
     # Fast-path rejection: if the array is empty or the target is out of range
     if isempty(mz_array) || target_mass + tolerance < first(mz_array) || target_mass - tolerance > last(mz_array)
@@ -930,59 +684,91 @@ function get_mz_slice(data::MSIData, mass::Real, tolerance::Real; mask_path::Uni
         precompute_analytics(data)
     end
 
-    println("Using high-performance sequential iterator...")
     target_min = mass - tolerance
     target_max = mass + tolerance
+    
+    # PERFORMANCE: Access raw vectors from metadata to avoid DataFrame dependency
+    # If stats_df exists, use it; otherwise, use the pre-computed bounds in SpectrumAsset
     stats_df = get_spectrum_stats(data)
+    
+    n_total = length(data.spectra_metadata)
+    # Instantiate thread-local buffer locally since global pools are deprecated
+    candidate_indices = Vector{Int}(undef, n_total)
+    
+    # We'll use local views of min/max if stats_df is missing to represent zero-allocation fallback
+    # But for extreme performance, we avoid list comprehensions [m... for m in ...] as they allocate.
+    local min_mzs::Vector{Float64}
+    local max_mzs::Vector{Float64}
+
+    if stats_df !== nothing
+        min_mzs = stats_df.MinMZ
+        max_mzs = stats_df.MaxMZ
+    else
+        # Fallback path: extract to local buffers or use metadata directly in loop
+        # For now, let's assume stats_df is usually populated by precompute_analytics.
+        # If not, we'll access it directly inside the filter loop.
+    end
+    candidate_count = 0
+    indices_to_check = masked_indices === nothing ? (1:n_total) : masked_indices
+    discretization_factor = 100.0
     bloom_filters = get_bloom_filters(data)
     
-    # 1. Find all candidate spectra first for efficient filtering
-    candidate_indices = Set{Int}()
-    indices_to_check = masked_indices === nothing ? (1:length(data.spectra_metadata)) : masked_indices
-    
     for i in indices_to_check
-        # NEW: Bloom filter check with discretization
-        if bloom_filters !== nothing && !is_empty(bloom_filters[i])
-            discretization_factor = 100.0
+        # Range check first (cheapest)
+        # Access metadata directly if stats_df is missing to ensure zero-allocation
+        @inbounds meta = data.spectra_metadata[i]
+        s_min, s_max = (stats_df !== nothing) ? (min_mzs[i], max_mzs[i]) : (meta.mz_asset.min_val, meta.mz_asset.max_val)
+        
+        if target_max < s_min || target_min > s_max
+            continue
+        end
+
+        # Bloom filter rejection (very fast)
+        if bloom_filters !== nothing
+            bf = bloom_filters[i]
+            if !is_empty(bf)
                 min_mass_int = round(Int, (mass - tolerance) * discretization_factor)
                 max_mass_int = round(Int, (mass + tolerance) * discretization_factor)
                 
                 found = false
-            for mass_int in min_mass_int:max_mass_int
-                if mass_int in bloom_filters[i]
+                @inbounds for mass_int in min_mass_int:max_mass_int
+                    if mass_int in bf
                         found = true
                         break
                     end
                 end
-            
-            if !found
-                continue # Definitely not in this spectrum
+                !found && continue
             end
         end
 
-        spec_min_mz = stats_df.MinMZ[i]
-        spec_max_mz = stats_df.MaxMZ[i]
-        if target_max >= spec_min_mz && target_min <= spec_max_mz
-            push!(candidate_indices, i)
-        end
+        candidate_count += 1
+        @inbounds candidate_indices[candidate_count] = i
     end
 
-    println("Found $(length(candidate_indices)) candidate spectra (filtered from $(length(indices_to_check)) initial spectra)")
+    # Use a view of the pre-allocated vector to avoid collect() allocations
+    valid_candidates = view(candidate_indices, 1:candidate_count)
     
     # 2. Iterate using the optimized, low-allocation iterator
-    results_count = 0
-    _iterate_spectra_fast(data, collect(candidate_indices)) do idx, mz_array, intensity_array
-        meta = data.spectra_metadata[idx]
+    # Use Atomic for thread-safe increment and avoid Ref-boxing
+    results_count = Base.Threads.Atomic{Int}(0)
+    
+    # Use let block to ensure closure captures are optimized (avoid boxing)
+    let slice_matrix=slice_matrix, results_count=results_count, width=width, height=height, 
+        spectra_metadata=data.spectra_metadata, mass=mass, tolerance=tolerance
+        
+        _iterate_spectra_fast(data, valid_candidates) do idx, mz_array, intensity_array
+            @inbounds meta = spectra_metadata[idx]
             intensity = find_mass(mz_array, intensity_array, mass, tolerance)
             if intensity > 0.0
                 if 1 <= meta.x <= width && 1 <= meta.y <= height
-                slice_matrix[meta.y, meta.x] = intensity
-                results_count += 1
+                    @inbounds slice_matrix[meta.y, meta.x] = intensity
+                    Base.Threads.atomic_add!(results_count, 1)
+                end
             end
         end
     end
     
-    println("Populated $results_count pixels with intensity data")
+    #println("Populated $(results_count[]) pixels with intensity data")
     replace!(slice_matrix, NaN => 0.0)
     return slice_matrix
 end
@@ -1007,13 +793,14 @@ This is a highly performant function that iterates through the full dataset only
 function get_multiple_mz_slices(data::MSIData, masses::AbstractVector{<:Real}, tolerance::Real; mask_path::Union{String, Nothing}=nothing)
     width, height = data.image_dims
     
-    # Sort masses to improve cache locality during search
+    # Sort masses to improve cache locality and allow binary search
     sorted_masses = sort(masses)
+    n_masses = length(sorted_masses)
 
-    # 1. Initialize a dictionary to hold the output slice matrices
-    slice_dict = Dict{Real, Matrix{Float64}}()
+    # 1. Initialize a dictionary to hold the output slice matrices (using Float32 for 50% RAM savings)
+    slice_dict = Dict{Real, Matrix{Float32}}()
     for mass in sorted_masses
-        slice_dict[mass] = zeros(Float64, height, width)
+        slice_dict[mass] = zeros(Float32, height, width)
     end
 
     local masked_indices::Union{Set{Int}, Nothing} = nothing
@@ -1029,42 +816,50 @@ function get_multiple_mz_slices(data::MSIData, masses::AbstractVector{<:Real}, t
         precompute_analytics(data)
     end
 
-    println("Filtering candidate spectra for $(length(masses)) m/z values...")
+    println("Filtering candidate spectra for $n_masses m/z values...")
     stats_df = get_spectrum_stats(data)
     bloom_filters = get_bloom_filters(data)
-    candidate_indices = Set{Int}()
+    
+    # Use a BitSet for faster index tracking
+    candidate_indices = BitSet()
     indices_to_check = masked_indices === nothing ? (1:length(data.spectra_metadata)) : masked_indices
 
-    # 3. Find all spectra that could contain *any* of the requested masses.
-    for mass in sorted_masses
-        target_min = mass - tolerance
-        target_max = mass + tolerance
-        for i in indices_to_check
-            # If already a candidate, no need to check again
-            if i in candidate_indices
-                continue
-            end
-            # NEW: Bloom filter check with discretization
-            if bloom_filters !== nothing && !is_empty(bloom_filters[i])
+    # 3. Optimized filtering: Iterate through spectra ONCE and check against all masses
+    # This changes complexity from O(M*N) to O(N * log M) or O(N + M) depending on range overlap
     discretization_factor = 100.0
+    for i in indices_to_check
+        spec_min = stats_df.MinMZ[i]
+        spec_max = stats_df.MaxMZ[i]
+        
+        # Binary search to find masses that might overlap with this spectrum's range
+        # target_min = mass - tolerance => mass = target_min + tolerance
+        # We need mass such that mass + tolerance >= spec_min  => mass >= spec_min - tolerance
+        # and mass - tolerance <= spec_max => mass <= spec_max + tolerance
+        
+        m_start_idx = searchsortedfirst(sorted_masses, spec_min - tolerance)
+        m_end_idx = searchsortedlast(sorted_masses, spec_max + tolerance)
+        
+        if m_start_idx <= m_end_idx
+            # Range overlap found, now check Bloom filter if available
+            if bloom_filters !== nothing && !is_empty(bloom_filters[i])
+                found_any = false
+                @inbounds for m_idx in m_start_idx:m_end_idx
+                    mass = sorted_masses[m_idx]
                     min_mass_int = round(Int, (mass - tolerance) * discretization_factor)
                     max_mass_int = round(Int, (mass + tolerance) * discretization_factor)
                     
-                found = false
                     for mass_int in min_mass_int:max_mass_int
                         if mass_int in bloom_filters[i]
-                        found = true
+                            found_any = true
                             break
                         end
                     end
-                
-                if !found
-                    continue # Definitely not in this spectrum
+                    found_any && break
                 end
+                if found_any
+                    push!(candidate_indices, i)
                 end
-            spec_min_mz = stats_df.MinMZ[i]
-            spec_max_mz = stats_df.MaxMZ[i]
-            if target_max >= spec_min_mz && target_min <= spec_max_mz
+            else
                 push!(candidate_indices, i)
             end
         end
@@ -1073,29 +868,40 @@ function get_multiple_mz_slices(data::MSIData, masses::AbstractVector{<:Real}, t
     println("Found $(length(candidate_indices)) total candidate spectra.")
 
     # 4. Iterate through the data a single time using the optimized iterator.
+    # We collect candidate_indices to pass to parallel iterator
     _iterate_spectra_fast(data, collect(candidate_indices)) do idx, mz_array, intensity_array
         meta = data.spectra_metadata[idx]
-        # For this single spectrum, check all masses of interest
-        for mass in sorted_masses
-            # Check if this spectrum's range actually covers the current mass
-            # This is a finer-grained check than the initial filtering
-            if !isempty(mz_array) && (mass + tolerance) >= first(mz_array) && (mass - tolerance) <= last(mz_array)
+        if isempty(mz_array)
+            return
+        end
+        
+        # Spectrum-level boundaries
+        spec_first = first(mz_array)
+        spec_last = last(mz_array)
+        
+        # Find which of our target masses fall within this specific spectrum's actual range
+        m_start_idx = searchsortedfirst(sorted_masses, spec_first - tolerance)
+        m_end_idx = searchsortedlast(sorted_masses, spec_last + tolerance)
+        
+        @inbounds for m_idx in m_start_idx:m_end_idx
+            mass = sorted_masses[m_idx]
             intensity = find_mass(mz_array, intensity_array, mass, tolerance)
             if intensity > 0.0
                 if 1 <= meta.x <= width && 1 <= meta.y <= height
-                        slice_dict[mass][meta.y, meta.x] = intensity
-                    end
+                    # Note: Concurrent writes to different matrices/coordinates are safe.
+                    # Dictionary access is safe because it's read-only after initialization.
+                    slice_dict[mass][meta.y, meta.x] = Float32(intensity)
                 end
             end
         end
     end
     
-    # 5. Clean up and return
+    # 5. Clean up - replaces NaNs with 0.0 directly in Float32 matrices
     for mass in sorted_masses
-        replace!(slice_dict[mass], NaN => 0.0)
+        replace!(slice_dict[mass], NaN32 => 0.0f0)
     end
     
-    println("Finished generating $(length(masses)) slices in a single pass.")
+    println("Finished generating $n_masses slices in a single pass.")
     return slice_dict
 end
 
@@ -1733,7 +1539,7 @@ Generates a colorbar image for a given slice of data.
 - `fig::Figure`: A figure with the colorbar.
 """
 function generate_colorbar_image(slice_data::AbstractMatrix, color_levels::Int, output_path::String, 
-                               bounds::Tuple{Float64, Float64}; 
+                               bounds::Tuple{Real, Real}; 
                                use_triq::Bool=false, triq_prob::Float64=0.98, 
                                mask_path::Union{String, Nothing}=nothing)
     # Use the provided bounds instead of recalculating

src/mzML.jl

@@ -231,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
@@ -394,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 ---")
-        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
+        seekstart(primary_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
@@ -433,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 ---
-        modes = [meta.mode for meta in spectra_metadata]
-        num_centroid = count(m -> m == CENTROID, modes)
-        num_profile = count(m -> m == PROFILE, modes)
+        # Determine overall acquisition mode ...
+        num_centroid = count(m -> m.mode == CENTROID, spectra_metadata)
+        num_profile = count(m -> m.mode == PROFILE, spectra_metadata)
         
         acq_mode_symbol = if num_centroid > 0 && num_profile == 0
             :centroid
@@ -468,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

@@ -19,6 +19,52 @@ 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_LTP/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_LA-ESI/180817_NEG_Thaliana_Leaf_bottom_1_0841.imzML",116.07,0.1, "Thaliana Leaf"),
+    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/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/imzML_AP_SMALDI/HR2MSImouseurinarybladderS096.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 MASK_ROUTE = "/home/pixel/Documents/Cinvestav_2025/JuliaMSI/public/css/masks/Stomach_DHB_uncompressed.png"
-# const MASK_ROUTE = "" 
+# const MASK_ROUTE = "/home/pixel/Documents/Cinvestav_2025/JuliaMSI/public/css/masks/Stomach_DHB_uncompressed.png"
+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_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()