.github/workflows/ci.yml

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

.gitignore

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

Dockerfile.headless

@@ -1,31 +0,0 @@
-# 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,5 +1,4 @@
 name = "MSI_src"
-uuid = "8f395f85-3b9c-4f7f-bf7e-12345678abcd"
 authors = ["JJSA"]
 version = "0.1.0"
 
@@ -36,7 +35,6 @@ 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"
@@ -44,11 +42,9 @@ 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]
@@ -56,9 +52,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"
@@ -67,6 +63,7 @@ GLMakie = "0.11"
 Genie = "5.31"
 GenieFramework = "2.8"
 HistogramThresholding = "0.3"
+Interpolations = "0.15"
 ImageBinarization = "0.3"
 ImageComponentAnalysis = "0.2"
 ImageContrastAdjustment = "0.3"
@@ -75,7 +72,6 @@ ImageFiltering = "0.7"
 ImageMorphology = "0.4"
 ImageSegmentation = "1.9"
 Images = "0.26"
-Interpolations = "0.15"
 JLD2 = "0.6"
 JSON = "0.21"
 Libz = "1.0"
@@ -94,12 +90,5 @@ 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,47 +34,6 @@ 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,9 +29,6 @@ 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()
@@ -63,7 +60,7 @@ if get(ENV, "GENIE_ENV", "dev") != "prod"
         println("--- MEMORY LOG [$(context)] ---")
         println("  Timestamp: $(now())")
         println("  Process RSS: $(rss_mb) MB")
-        println("  global_msi_data[] size: $(msi_data_size_mb) MB")
+        println("  msi_data size: $(msi_data_size_mb) MB")
         println("  Cumulative GC time: $(gc_time_s) s")
         println("--------------------------")
     end
@@ -196,15 +193,6 @@ 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 ==
@@ -488,7 +476,7 @@ Core.eval(Stipple, :(sesstoken() = "")) # Prevent ErrorException("Model storage
 
     # == DATA MANAGEMENT VARIABLES ==
     # Centralized MSIData object
-    # global_msi_data[] is now global to avoid Genie Session memory leaks
+    @out msi_data::Union{MSIData, Nothing} = nothing
 
     # Image file management
     @out text_nmass="" # For specific mass charge image creation
@@ -673,7 +661,7 @@ Core.eval(Stipple, :(sesstoken() = "")) # Prevent ErrorException("Model storage
         
         try
             # 1. Clear large data objects explicitly
-            global_msi_data[] = nothing
+            msi_data = nothing
             feature_matrix_result = nothing
             bin_info_result = nothing
             
@@ -729,42 +717,30 @@ Core.eval(Stipple, :(sesstoken() = "")) # Prevent ErrorException("Model storage
     sure the file can be processed by later steps like mainProcess 
     =#
     @onbutton btnSearch begin
-        # 0. Robustness Guard: Prevent double-trigger during processing
+        is_processing = true
-        if is_processing
+        push!(__model__)
-            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
+        # --- Close previous dataset if one is open ---
-        is_processing = true
+        if msi_data !== nothing
+            println("DEBUG: Closing previously loaded dataset before opening new one: $(basename(full_route))")
+            close(msi_data)
+            msi_data = nothing
+            GC.gc()
+            if Sys.islinux()
+                ccall(:malloc_trim, Int32, (Int32,), 0)
+            end
+        end
+
         msg = "Opening file: $(basename(picked_route))..."
-        SpectraEnabled = false
-        btnMetadataDisable = true
-        push!(__model__)
 
-        # 3. Spawn background computational thread
+        try
-        Threads.@spawn begin
+            dataset_name = replace(basename(picked_route), r"(\.(imzML|imzml|mzML|mzml))$"i => "")
-            try
-                # --- Close previous dataset if one is open ---
-                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
-
-                dataset_name = replace(basename(picked_route), r"(\.(imzML|imzml|mzML|mzml))$"i => "")
             registry = load_registry(registry_path)
             existing_entry = get(registry, dataset_name, nothing)
 
@@ -781,8 +757,8 @@ Core.eval(Stipple, :(sesstoken() = "")) # Prevent ErrorException("Model storage
                 dims = parse.(Int, split(dims_str, " x "))
                 imgWidth, imgHeight = dims[1], dims[2]
 
-                global_msi_data[] = nothing # Ensure data is not held in memory
+                msi_data = nothing # Ensure data is not held in memory
-                log_memory_usage("Fast Load (global_msi_data[] cleared)", global_msi_data[])
+                log_memory_usage("Fast Load (msi_data cleared)", msi_data)
                 btnMetadataDisable = false
                 SpectraEnabled = true
                 selected_folder_main = dataset_name
@@ -1024,10 +1000,10 @@ Core.eval(Stipple, :(sesstoken() = "")) # Prevent ErrorException("Model storage
             image_available_folders = deepcopy(img_folders)
 
             selected_folder_main = dataset_name
-            global_msi_data[] = loaded_data
+            msi_data = loaded_data
 
             # Determine plot mode from loaded data
-            df = global_msi_data[].spectrum_stats_df
+            df = 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)
@@ -1037,29 +1013,26 @@ Core.eval(Stipple, :(sesstoken() = "")) # Prevent ErrorException("Model storage
                 last_plot_mode = "lines" # Default
             end
 
-            log_memory_usage("Full Load", global_msi_data[])
+            log_memory_usage("Full Load", msi_data)
 
             eTime = round(time() - sTime, digits=3)
             msg = "Active file loaded in $(eTime) seconds. Dataset '$(dataset_name)' is ready for analysis."
 
             SpectraEnabled = true
             
-            catch e
+        catch e
-                global_msi_data[] = nothing
+            msi_data = nothing
-                msg = "Error loading active file: $e"
+            msg = "Error loading active file: $e"
-                warning_msg = true
+            warning_msg = true
-                SpectraEnabled = false
+            SpectraEnabled = false
-                btnMetadataDisable = true
+            btnMetadataDisable = true
-                @error "File loading failed" exception=(e, catch_backtrace())
+            @error "File loading failed" exception=(e, catch_backtrace())
-                push!(__model__) # Force sending error back to UI immediately
+        finally
-            finally
+            GC.gc() # Trigger garbage collection
-                GC.gc() # Trigger garbage collection
+            if Sys.islinux()
-                if Sys.islinux()
+                ccall(:malloc_trim, Int32, (Int32,), 0) # Ensure Julia returns the freed memory to OS
-                    ccall(:malloc_trim, Int32, (Int32,), 0) # Ensure Julia returns the freed memory to OS
-                end
-                is_processing = false
-                push!(__model__) # Clear spinner loop 
             end
+            is_processing = false
         end
     end
 
@@ -1259,11 +1232,6 @@ Core.eval(Stipple, :(sesstoken() = "")) # Prevent ErrorException("Model storage
     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
@@ -1291,9 +1259,9 @@ Core.eval(Stipple, :(sesstoken() = "")) # Prevent ErrorException("Model storage
             end
             target_path = entry["source_path"]
 
-            # Ensure global_msi_data[] is for the currently selected file and load if needed
+            # Ensure 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 global_msi_data[] used for plotting/interactive exploration.
+            # with the global msi_data used for plotting/interactive exploration.
             println("DEBUG: Opening isolated MSIData instance for pipeline stability...")
             pipeline_msi_data = OpenMSIData(target_path)
 
@@ -1718,9 +1686,9 @@ Core.eval(Stipple, :(sesstoken() = "")) # Prevent ErrorException("Model storage
                     
                     # Determine plot mode for this specific spectrum
                     spectrum_mode_for_plot = "lines" # Default to lines
-                    if global_msi_data[].spectrum_stats_df !== nothing && "Mode" in names(global_msi_data[].spectrum_stats_df)
+                    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(global_msi_data[].spectrum_stats_df.Mode)
+                        if selected_spectrum_id_for_plot > 0 && selected_spectrum_id_for_plot <= length(msi_data.spectrum_stats_df.Mode)
-                            mode = global_msi_data[].spectrum_stats_df.Mode[selected_spectrum_id_for_plot]
+                            mode = msi_data.spectrum_stats_df.Mode[selected_spectrum_id_for_plot]
                             if mode == MSI_src.CENTROID
                                 spectrum_mode_for_plot = "stem"
                             end
@@ -1971,7 +1939,7 @@ Core.eval(Stipple, :(sesstoken() = "")) # Prevent ErrorException("Model storage
     @onbutton recalculate_suggestions_btn begin
         is_processing = true
         push!(__model__)
-        if global_msi_data[] === nothing
+        if msi_data === nothing
             msg = "Please load a file first."
             warning_msg = true
             return
@@ -1985,7 +1953,7 @@ Core.eval(Stipple, :(sesstoken() = "")) # Prevent ErrorException("Model storage
                 for p in reference_peaks_list 
                 if tryparse(Float64, string(p["mz"])) !== nothing
             )
-            recommended_params = main_precalculation(global_msi_data[], reference_peaks=ref_peaks)
+            recommended_params = main_precalculation(msi_data, reference_peaks=ref_peaks)
             
 
             for (step_name, params) in recommended_params
@@ -2380,11 +2348,6 @@ Core.eval(Stipple, :(sesstoken() = "")) # Prevent ErrorException("Model storage
     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..."
@@ -2459,22 +2422,8 @@ Core.eval(Stipple, :(sesstoken() = "")) # Prevent ErrorException("Model storage
                 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)
@@ -2609,21 +2558,21 @@ Core.eval(Stipple, :(sesstoken() = "")) # Prevent ErrorException("Model storage
                 return
             end
 
-            if global_msi_data[] === nothing || full_route != target_path
+            if msi_data === nothing || full_route != target_path
-                if global_msi_data[] !== nothing
+                if msi_data !== nothing
-                    close(global_msi_data[])
+                    close(msi_data)
                 end
                 msg = "Reloading $(basename(target_path)) for analysis..."
                 full_route = target_path
-                global_msi_data[] = OpenMSIData(target_path)
+                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!(global_msi_data[], convert(Float64, min_val), convert(Float64, max_val))
+                    set_global_mz_range!(msi_data, convert(Float64, min_val), convert(Float64, max_val))
                 else
-                    precompute_analytics(global_msi_data[])
+                    precompute_analytics(msi_data)
                 end
             end
 
@@ -2636,7 +2585,7 @@ Core.eval(Stipple, :(sesstoken() = "")) # Prevent ErrorException("Model storage
                 end
             end
 
-            plotdata, plotlayout, xSpectraMz, ySpectraMz = meanSpectrumPlot(global_msi_data[], selected_folder_main, mask_path=mask_path_for_plot)
+            plotdata, plotlayout, xSpectraMz, ySpectraMz = meanSpectrumPlot(msi_data, selected_folder_main, mask_path=mask_path_for_plot)
             plotdata_before = plotdata
             plotlayout_before = plotlayout
             last_plot_type = "mean"
@@ -2644,7 +2593,7 @@ Core.eval(Stipple, :(sesstoken() = "")) # Prevent ErrorException("Model storage
             fTime = time()
             eTime = round(fTime - sTime, digits=3)
             msg = "Plot loaded in $(eTime) seconds"
-            log_memory_usage("Mean Plot Generated", global_msi_data[])
+            log_memory_usage("Mean Plot Generated", msi_data)
         catch e
             msg = "Could not generate mean spectrum plot: $e"
             warning_msg = true
@@ -2698,21 +2647,21 @@ Core.eval(Stipple, :(sesstoken() = "")) # Prevent ErrorException("Model storage
                 return
             end
 
-            if global_msi_data[] === nothing || full_route != target_path
+            if msi_data === nothing || full_route != target_path
-                if global_msi_data[] !== nothing
+                if msi_data !== nothing
-                    close(global_msi_data[])
+                    close(msi_data)
                 end
                 msg = "Reloading $(basename(target_path)) for analysis..."
                 full_route = target_path
-                global_msi_data[] = OpenMSIData(target_path)
+                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!(global_msi_data[], convert(Float64, min_val), convert(Float64, max_val))
+                    set_global_mz_range!(msi_data, convert(Float64, min_val), convert(Float64, max_val))
                 else
-                    precompute_analytics(global_msi_data[])
+                    precompute_analytics(msi_data)
                 end
             end
             local mask_path_for_plot::Union{String, Nothing} = nothing
@@ -2724,7 +2673,7 @@ Core.eval(Stipple, :(sesstoken() = "")) # Prevent ErrorException("Model storage
                 end
             end
 
-            plotdata, plotlayout, xSpectraMz, ySpectraMz = sumSpectrumPlot(global_msi_data[], selected_folder_main, mask_path=mask_path_for_plot)
+            plotdata, plotlayout, xSpectraMz, ySpectraMz = sumSpectrumPlot(msi_data, selected_folder_main, mask_path=mask_path_for_plot)
             plotdata_before = plotdata
             plotlayout_before = plotlayout
             last_plot_type = "sum"
@@ -2732,7 +2681,7 @@ Core.eval(Stipple, :(sesstoken() = "")) # Prevent ErrorException("Model storage
             fTime = time()
             eTime = round(fTime - sTime, digits=3)
             msg = "Total plot loaded in $(eTime) seconds"
-            log_memory_usage("Sum Plot Generated", global_msi_data[])
+            log_memory_usage("Sum Plot Generated", msi_data)
         catch e
             msg = "Could not generate total spectrum plot: $e"
             warning_msg = true
@@ -2789,21 +2738,21 @@ Core.eval(Stipple, :(sesstoken() = "")) # Prevent ErrorException("Model storage
                 return
             end
 
-            if global_msi_data[] === nothing || full_route != target_path
+            if msi_data === nothing || full_route != target_path
-                if global_msi_data[] !== nothing
+                if msi_data !== nothing
-                    close(global_msi_data[])
+                    close(msi_data)
                 end
                 msg = "Reloading $(basename(target_path)) for analysis..."
                 full_route = target_path
-                global_msi_data[] = OpenMSIData(target_path)
+                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!(global_msi_data[], convert(Float64, min_val), convert(Float64, max_val))
+                    set_global_mz_range!(msi_data, convert(Float64, min_val), convert(Float64, max_val))
                 else
-                    precompute_analytics(global_msi_data[])
+                    precompute_analytics(msi_data)
                 end
             end
 
@@ -2818,7 +2767,7 @@ Core.eval(Stipple, :(sesstoken() = "")) # Prevent ErrorException("Model storage
 
             # Convert to positive coordinates for processing
             y_positive = yCoord < 0 ? abs(yCoord) : yCoord
-            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, plotlayout, xSpectraMz, ySpectraMz, spectrum_id = xySpectrumPlot(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"
@@ -2859,7 +2808,7 @@ Core.eval(Stipple, :(sesstoken() = "")) # Prevent ErrorException("Model storage
             fTime = time()
             eTime = round(fTime - sTime, digits=3)
             msg = "Plot loaded in $(eTime) seconds"
-            log_memory_usage("XY Plot Generated", global_msi_data[])
+            log_memory_usage("XY Plot Generated", msi_data)
         catch e
             msg = "Could not retrieve spectrum: $e"
             warning_msg = true
@@ -2904,21 +2853,21 @@ Core.eval(Stipple, :(sesstoken() = "")) # Prevent ErrorException("Model storage
                 return
             end
 
-            if global_msi_data[] === nothing || full_route != target_path
+            if msi_data === nothing || full_route != target_path
-                if global_msi_data[] !== nothing
+                if msi_data !== nothing
-                    close(global_msi_data[])
+                    close(msi_data)
                 end
                 msg = "Reloading $(basename(target_path)) for analysis..."
                 full_route = target_path
-                global_msi_data[] = OpenMSIData(target_path)
+                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!(global_msi_data[], convert(Float64, min_val), convert(Float64, max_val))
+                    set_global_mz_range!(msi_data, convert(Float64, min_val), convert(Float64, max_val))
                 else
-                    precompute_analytics(global_msi_data[])
+                    precompute_analytics(msi_data)
                 end
             end
 
@@ -2932,7 +2881,7 @@ Core.eval(Stipple, :(sesstoken() = "")) # Prevent ErrorException("Model storage
             end
 
             # Call the new nSpectrumPlot function
-            plotdata, plotlayout, xSpectraMz, ySpectraMz, spectrum_id = nSpectrumPlot(global_msi_data[], idSpectrum, selected_folder_main, mask_path=mask_path_for_plot)
+            plotdata, plotlayout, xSpectraMz, ySpectraMz, spectrum_id = nSpectrumPlot(msi_data, idSpectrum, selected_folder_main, mask_path=mask_path_for_plot)
             plotdata_before = plotdata
             plotlayout_before = plotlayout
             last_plot_type = "single"
@@ -2942,7 +2891,7 @@ Core.eval(Stipple, :(sesstoken() = "")) # Prevent ErrorException("Model storage
             fTime = time()
             eTime = round(fTime - sTime, digits=3)
             msg = "Plot loaded in $(eTime) seconds"
-            log_memory_usage("nSpectrum Plot Generated", global_msi_data[])
+            log_memory_usage("nSpectrum Plot Generated", msi_data)
         catch e
             msg = "Could not retrieve spectrum: $e"
             warning_msg = true
@@ -3223,7 +3172,7 @@ Core.eval(Stipple, :(sesstoken() = "")) # Prevent ErrorException("Model storage
 
     # This handler will now correctly load the first image from the newly selected folder.
     @onchange selected_folder_main begin
-        # The global_msi_data[] object lifecycle is managed by the btnSearch handler.
+        # The 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)
@@ -3455,7 +3404,7 @@ Core.eval(Stipple, :(sesstoken() = "")) # Prevent ErrorException("Model storage
             fTime = time()
             eTime = round(fTime - sTime, digits=3)
             msg = "Plot loaded in $(eTime) seconds"
-            log_memory_usage("Mean Plot Generated", global_msi_data[])
+            log_memory_usage("Mean Plot Generated", msi_data)
         catch e
             msg = "Failed to load and process image: $e"
             warning_msg = true
@@ -3513,7 +3462,7 @@ Core.eval(Stipple, :(sesstoken() = "")) # Prevent ErrorException("Model storage
             fTime = time()
             eTime = round(fTime - sTime, digits=3)
             msg = "Plot loaded in $(eTime) seconds"
-            log_memory_usage("Mean Plot Generated", global_msi_data[])
+            log_memory_usage("Mean Plot Generated", msi_data)
         catch e
             msg = "Failed to load and process image: $e"
             warning_msg = true
@@ -3614,7 +3563,7 @@ Core.eval(Stipple, :(sesstoken() = "")) # Prevent ErrorException("Model storage
     # To include a visualization in the spectrum plot indicating where is the selected mass
     @onchange Nmass begin
         if !isempty(xSpectraMz)
-            df = global_msi_data[].spectrum_stats_df
+            df = 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)
@@ -3873,7 +3822,7 @@ Core.eval(Stipple, :(sesstoken() = "")) # Prevent ErrorException("Model storage
         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", global_msi_data[])
+        log_memory_usage("App Ready", msi_data)
     end
     # is_processing = false
     GC.gc() # Trigger garbage collection

app.jl.html

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

build_sysimage.jl

@@ -1,49 +0,0 @@
-# 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,18 +903,6 @@ 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)
 
@@ -926,9 +914,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::String)
+function load_registry(registry_path)
     return lock(REGISTRY_LOCK) do
-        clean_path = clean_registry_path(registry_path)
+        clean_path = abspath(registry_path)
         if isfile(clean_path)
             for attempt in 1:5
                 try
@@ -949,82 +937,93 @@ function load_registry(registry_path::String)
 end
 
 """
-    extract_metadata(loaded_data::MSIData, full_route::String)
+    extract_metadata(msi_data::MSIData, source_path::String)
 
-Extracts key metadata from an MSIData object and file path for UI display.
+Extracts key metadata from an `MSIData` object for display.
-Returns a dictionary with a "summary" key containing parameter-value pairs.
+
+# Arguments
+- `msi_data`: The `MSIData` object.
+- `source_path`: The path to the source data file.
+
+# Returns
+- A dictionary containing summary statistics and metadata.
 """
-function extract_metadata(loaded_data::MSIData, full_route::String)
+function extract_metadata(msi_data::MSIData, source_path::String)
-    stats = []
+    df = msi_data.spectrum_stats_df
-    
+    if df === nothing
-    # Basic File Info
+        # This can happen if precompute_analytics hasn't been run
-    push!(stats, Dict("parameter" => "Filename", "value" => basename(full_route)))
+        # We can still return basic info
-    push!(stats, Dict("parameter" => "Source Path", "value" => full_route))
+        return Dict(
-    
+            "summary" => [
-    # Image Dimensions
+                Dict("parameter" => "File Name", "value" => basename(source_path)),
-    w, h = loaded_data.image_dims
+                Dict("parameter" => "Number of Spectra", "value" => length(msi_data.spectra_metadata)),
-    if w > 0 && h > 0
+                Dict("parameter" => "Image Dimensions", "value" => "$(msi_data.image_dims[1]) x $(msi_data.image_dims[2])"),
-        push!(stats, Dict("parameter" => "Image Dimensions", "value" => "$w x $h"))
+            ],
-    else
+            "global_min_mz" => nothing,
-        push!(stats, Dict("parameter" => "Image Dimensions", "value" => "N/A (Non-imaging)"))
+            "global_max_mz" => nothing
+        )
     end
-    
+
-    # Spectrum Counts
+    summary_stats = [
-    num_spectra = length(loaded_data.spectra_metadata)
+        Dict("parameter" => "File Name", "value" => basename(source_path)),
-    push!(stats, Dict("parameter" => "Total Spectra", "value" => string(num_spectra)))
+        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])"),
-    # m/z Range
+        Dict("parameter" => "Global Min m/z", "value" => @sprintf("%.4f", Threads.atomic_add!(msi_data.global_min_mz, 0.0))),
-    min_mz, max_mz = get_global_mz_range(loaded_data)
+        Dict("parameter" => "Global Max m/z", "value" => @sprintf("%.4f", Threads.atomic_add!(msi_data.global_max_mz, 0.0))),
-    if isfinite(min_mz) && isfinite(max_mz)
+        Dict("parameter" => "Mean TIC", "value" => @sprintf("%.2e", mean(df.TIC))),
-        push!(stats, Dict("parameter" => "Global m/z Range", "value" => "$(round(min_mz, digits=4)) - $(round(max_mz, digits=4))"))
+        Dict("parameter" => "Mean BPI", "value" => @sprintf("%.2e", mean(df.BPI))),
-    else
+        Dict("parameter" => "Mean # Points", "value" => @sprintf("%.1f", mean(df.NumPoints))),
-        push!(stats, Dict("parameter" => "Global m/z Range", "value" => "Unknown"))
+    ]
-    end
+
-    
+    if hasproperty(df, :Mode)
-    # Instrument Metadata
+        centroid_count = count(==(MSI_src.CENTROID), df.Mode)
-    if loaded_data.instrument_metadata !== nothing
+        profile_count = count(==(MSI_src.PROFILE), df.Mode)
-        inst = loaded_data.instrument_metadata
+        unknown_count = count(==(MSI_src.UNKNOWN), df.Mode)
-        push!(stats, Dict("parameter" => "Instrument Model", "value" => inst.instrument_model))
+
-        push!(stats, Dict("parameter" => "Polarity", "value" => titlecase(string(inst.polarity))))
+        push!(summary_stats, Dict("parameter" => "Centroid Spectra", "value" => string(centroid_count)))
-        push!(stats, Dict("parameter" => "Acquisition Mode", "value" => titlecase(string(inst.acquisition_mode))))
+        push!(summary_stats, Dict("parameter" => "Profile Spectra", "value" => string(profile_count)))
-        
+        if unknown_count > 0
-        if inst.resolution !== nothing
+            push!(summary_stats, Dict("parameter" => "Unknown Mode Spectra", "value" => string(unknown_count)))
-            push!(stats, Dict("parameter" => "Instrument Resolution", "value" => string(inst.resolution)))
         end
     end
-    
+
-    return Dict("summary" => stats)
+    return Dict(
+        "summary" => summary_stats,
+        "global_min_mz" => msi_data.global_min_mz,
+        "global_max_mz" => msi_data.global_max_mz
+    )
 end
 
 """
-    save_registry(registry_path::String, registry_data::Dict)
+    save_registry(registry_path, registry_data)
 
-Saves the dataset registry to a JSON file atomically to prevent data loss 
+Saves the dataset registry to a JSON file, ensuring thread-safe access and robustness on Windows.
-and sharing violations on Windows. Writes to a `.tmp` file first.
+
+# Arguments
+- `registry_path`: Path to the `registry.json` file.
+- `registry_data`: The dictionary containing the registry data to save.
 """
-function save_registry(registry_path::String, registry_data)
+function save_registry(registry_path, registry_data)
     lock(REGISTRY_LOCK) do
-        clean_path = clean_registry_path(registry_path)
+        # Ensure the path is clean and valid for the OS
+        clean_path = abspath(registry_path)
         dir = dirname(clean_path)
         if !isdir(dir)
             mkpath(dir)
         end
 
-        temp_path = clean_path * ".tmp"
+        # On Windows, opening files for writing can occasionally fail if another process 
-
+        # (like an anti-virus or file indexer) holds a transient lock.
-        for attempt in 1:6
+        # Retry with a small delay to improve robustness.
+        for attempt in 1:5
             try
-                open(temp_path, "w") do f
+                open(clean_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 == 6
+                if attempt == 5
                     @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..."
@@ -1048,21 +1047,29 @@ 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::String, dataset_name::String, source_path::String, metadata=nothing, is_imzML=false)
+function update_registry(registry_path, dataset_name, source_path, 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
@@ -1131,24 +1138,23 @@ function process_file_safely(file_path, masses, params, progress_message_ref, ov
         # --- Extract metadata ---
         metadata = extract_metadata(local_msi_data, file_path)
 
-        # --- Save Slices (Parallelized) ---
+        # --- Save Slices ---
         mkpath(output_dir)
-        
+
-        progress_lock = ReentrantLock()
+        for (mass_idx, mass) in enumerate(masses)
-        
+            progress_message_ref = "File $(params.fileIdx)/$(params.nFiles): Saving slice for m/z=$mass"
-        # Parallelize over the requested masses to fully utilize multi-core CPUs
+
-        Threads.@threads for mass in masses
             slice = slice_dict[mass]
             text_nmass = replace(string(mass), "." => "_")
             bitmap_filename = params.triqE ? "TrIQ_$(text_nmass).bmp" : "MSI_$(text_nmass).bmp"
             colorbar_filename = params.triqE ? "colorbar_TrIQ_$(text_nmass).png" : "colorbar_MSI_$(text_nmass).png"
 
-            local sliceQuant, bounds
             if all(iszero, slice)
                 sliceQuant = zeros(UInt8, size(slice))
                 bounds = (0.0, 1.0)  # Default bounds for empty slice
+                @warn "No intensity data for m/z = $mass in $(dataset_name)"
             else
-                # TrIQ and quantization are computationally intensive and now run in parallel
+                # Get both quantized data AND bounds in one call
                 if params.triqE
                     sliceQuant, bounds = TrIQ(slice, params.colorL, params.triqP, mask_matrix=mask_matrix_for_triq)
                 else
@@ -1157,24 +1163,18 @@ 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)
-                # CairoMakie colorbar generation is now also parallelized
+                # Now pass the precomputed bounds to colorbar generation
                 generate_colorbar_image(slice, params.colorL, joinpath(output_dir, colorbar_filename), 
                                       bounds; use_triq=params.triqE, triq_prob=params.triqP, mask_path=mask_path)
             end
-            
-            lock(progress_lock) do
-                progress_message_ref = "File $(params.fileIdx)/$(params.nFiles): Saved slice for m/z=$mass"
-            end
         end
 
-
         is_imzML = local_msi_data.source isa ImzMLSource
         update_registry(params.registry, dataset_name, file_path, metadata, is_imzML)
         return (true, "")
@@ -1194,14 +1194,24 @@ function process_file_safely(file_path, masses, params, progress_message_ref, ov
     end
 end
 
-function update_registry_mask_fields(registry_path::String, dataset_name::String, has_mask::Bool, mask_path::String)
+function update_registry_mask_fields(registry_path, dataset_name, has_mask, mask_path)
     lock(REGISTRY_LOCK) do
-        registry = load_registry(registry_path)
+        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
         
         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,
@@ -1211,6 +1221,7 @@ function update_registry_mask_fields(registry_path::String, dataset_name::String
             )
         end
 
+        # Use the robust save function
         save_registry(registry_path, registry)
     end
 end

mask.jl

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

precompile_script.jl

@@ -1,173 +0,0 @@
-# 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,7 +20,30 @@ 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,37 +1,77 @@
+# 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
 
+# --- Buffer Pooling ---
 """
-    posix_madvise(buffer::AbstractArray, advice::Integer)
+    BufferPool
 
-A safe wrapper for the OS `madvise` system call. Signals the kernel about the 
+A thread-safe pool of `Vector{UInt8}` buffers, categorized by their size.
-access pattern for a memory-mapped region. Currently supports Linux/Unix systems.
+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.
 """
-function posix_madvise(buffer::AbstractArray, advice::Integer)
+mutable struct BufferPool
-    @static if Sys.isunix()
+    lock::ReentrantLock
-        try
+    buffers::Dict{Int, Vector{Vector{UInt8}}}
-            ptr = pointer(buffer)
+
-            len = sizeof(buffer)
+    function BufferPool()
-            # ccall(:madvise, return_type, (arg_types...), args...)
+        new(ReentrantLock(), Dict{Int, Vector{Vector{UInt8}}}())
-            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 ---
+"""
+    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
 
 
 """
@@ -71,16 +111,11 @@ If no suitable buffer is available, a new one is allocated.
 - A `Vector{UInt8}` buffer.
 """
 function get_buffer!(pool::SimpleBufferPool, size::Int)::Vector{UInt8}
-    buf = lock(pool.lock) do
+    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,18 +20,7 @@ export OpenMSIData,
        MSIData,
        _iterate_spectra_fast,
        validate_spectrum,
-       get_mz_slice,
+       REGISTRY_LOCK
-       REGISTRY_LOCK,
-       SpectrumMetadata,
-       SpectrumAsset,
-       SpectrumMode,
-       CENTROID,
-       PROFILE,
-       MzMLSource,
-       ImzMLSource,
-       SimpleBufferPool,
-       get_buffer!,
-       release_buffer!
 
 # Define shared registry lock
 const REGISTRY_LOCK = ReentrantLock()
@@ -71,21 +60,9 @@ 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")
@@ -95,8 +72,6 @@ include("Preprocessing.jl")
 include("ImageProcessing.jl")
 include("Precalculations.jl")
 include("PreprocessingPipeline.jl")
-include("StreamingKernels.jl")
-include("StreamingPipeline.jl")
 
 using Setfield # For immutable struct updates
 
@@ -120,17 +95,13 @@ 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(norm_filepath), ".mzml")
+    if endswith(lowercase(filepath), ".mzml")
-        msi_data = load_mzml_lazy(norm_filepath, cache_size=cache_size)
+        msi_data = load_mzml_lazy(filepath, cache_size=cache_size)
-    elseif endswith(lowercase(norm_filepath), ".imzml")
+    elseif endswith(lowercase(filepath), ".imzml")
-        msi_data = load_imzml_lazy(norm_filepath, cache_size=cache_size)
+        msi_data = load_imzml_lazy(filepath, cache_size=cache_size)
     else
-        error("Unsupported file type: $norm_filepath. Please provide a .mzML or .imzML file.")
+        error("Unsupported file type: $filepath. Please provide a .mzML or .imzML file.")
     end
 
     # Apply spectrum type map if provided

src/Precalculations.jl

@@ -1314,16 +1314,18 @@ 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
-    count = end_idx - start_idx + 1
+    if (end_idx - start_idx + 1) < 3
-    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 = float(intensity[peak_idx])
+    A_est = intensity[peak_idx]
     # Mean (μ): m/z at peak intensity
-    mu_est = float(mz[peak_idx])
+    mu_est = 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)
@@ -1337,27 +1339,19 @@ function _fit_gaussian_and_r2(mz::AbstractVector{<:Real}, intensity::AbstractVec
         return 0.0
     end
 
-    # Calculate SS_res, mean_y in a single pass to avoid allocations
+    # Gaussian function
-    SS_res = 0.0
+    gaussian(x, A, mu, sigma) = A * exp.(-(x .- mu).^2 ./ (2 * sigma^2))
-    sum_y = 0.0
-    
-    @inbounds for i in start_idx:end_idx
-        x_val = float(mz[i])
-        y_val = float(intensity[i])
-        
-        # Gaussian function estimate
-        y_est = A_est * exp(-((x_val - mu_est)^2) / (2 * sigma_est^2))
-        
-        SS_res += (y_val - y_est)^2
-        sum_y += y_val
-    end
 
-    mean_y = sum_y / count
+    # Generate estimated Gaussian curve
-    SS_tot = 0.0
+    y_est = gaussian(x_data, A_est, mu_est, sigma_est)
-    
+
-    @inbounds for i in start_idx:end_idx
+    # Calculate pseudo R-squared
-        SS_tot += (float(intensity[i]) - mean_y)^2
+    # R^2 = 1 - (SS_res / SS_tot)
-    end
+    # SS_res = sum((y_data - y_est).^2)
+    # SS_tot = sum((y_data - mean(y_data)).^2)
+
+    SS_res = sum((y_data .- y_est).^2)
+    SS_tot = sum((y_data .- mean(y_data)).^2)
 
     if SS_tot == 0
         return 1.0 # Perfect fit if all y_data are the same

src/Preprocessing.jl

@@ -12,7 +12,6 @@ generation.
 # =============================================================================
 
 using Statistics  # For mean, median
-using SparseArrays
 using StatsBase   # For mad (Median Absolute Deviation)
 using SavitzkyGolay # For SavitzkyGolay filtering
 using Dates         # For now()
@@ -41,7 +40,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::AbstractMatrix{Float64}
+    matrix::Array{Float64,2}
     mz_bins::Vector{Tuple{Float64,Float64}}
     sample_ids::Vector{Int}
 end
@@ -488,23 +487,32 @@ 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 the baseline estimate array once
+    # 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
     b1 = collect(float.(y))
+    b2 = similar(b1)        
     
+    current_b = b1
+    next_b = b2
+
     for k in 1:iterations
-        prev_val = b1[1]
+        # Calculate next baseline estimate into `next_b` based on `current_b`
-        b1[1] = min(b1[1], b1[2])
+        # Boundary conditions
-        
+        if n > 1
-        @inbounds for i in 2:n-1
+            next_b[1] = min(current_b[1], current_b[2])
-            curr_val = b1[i]
+            next_b[n] = min(current_b[n], current_b[n-1])
-            b1[i] = min(curr_val, 0.5 * (prev_val + b1[i+1]))
-            prev_val = curr_val
         end
-        b1[n] = min(b1[n], prev_val)
+
+        @inbounds for i in 2:n-1
+            next_b[i] = min(current_b[i], 0.5 * (current_b[i-1] + current_b[i+1]))
+        end
+        
+        # Swap references for the next iteration (no data copy here)
+        current_b, next_b = next_b, current_b 
     end
     
-    # Return the final baseline estimate
+    # Return the final baseline estimate (which is in current_b after the last swap)
-    return b1
+    return current_b
 end
 
 """
@@ -712,32 +720,18 @@ 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}}[]
     
-    # Fast, non-allocating noise estimation
+    noise_level = mad(y, normalize=true) + eps(Float64)
-    mean_y = sum(y) / n
+    ys = smooth_spectrum_core(y; method=:savitzky_golay, window=max(5, 2*half_window+1), order=2) # Use smoothed data for detection
-    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[]
-    sizehint!(candidate_peak_indices, div(n, 10)) # Pre-allocate memory capacity
+    for i in 2:n-1
-    
-    @inbounds for i in 2:n-1
         left = max(1, i - half_window)
         right = min(n, i + half_window)
         
-        # Avoid @view allocation in tight loop by manually computing minimums and maximums
+        # Prominence check
-        min_left = ys[left]
+        prominence = ys[i] - max(minimum(@view ys[left:i]), minimum(@view ys[i:right]))
-        for j in left:i; min_left = min(min_left, ys[j]); end
         
-        min_right = ys[i]
+        if ys[i] >= maximum(@view ys[left:right]) && 
-        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)
@@ -774,14 +768,7 @@ 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

@@ -1,79 +0,0 @@
-# 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

@@ -1,317 +0,0 @@
-# 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

@@ -1,402 +0,0 @@
-# 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,4 +1,5 @@
-using Images, Statistics, CairoMakie, DataFrames, Printf, ColorSchemes, StatsBase, Mmap
+# src/imzML.jl
+using Images, Statistics, CairoMakie, DataFrames, Printf, ColorSchemes, StatsBase
 
 """
 This file provides a library for parsing `.imzML` and `.ibd` files in pure Julia.
@@ -194,6 +195,86 @@ 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)
 
@@ -374,7 +455,7 @@ end
 # Arguments
 
 - `stream::IO`: The stream to parse.
-- `hIbd::IO`: The IBD file handle.
+- `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.
@@ -389,7 +470,7 @@ end
 
 - `spectra_metadata::Vector{SpectrumMetadata}`: The spectrum metadata.
 """
-function parse_imzml_spectrum_block(stream::IO, hIbd::IO, param_groups::Dict{String, SpecDim},
+function parse_imzml_spectrum_block(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)
@@ -406,8 +487,8 @@ function parse_imzml_spectrum_block(stream::IO, hIbd::IO, param_groups::Dict{Str
             @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, 0.0, 0.0)
+                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, 0.0, 0.0)
+                int_asset = SpectrumAsset(default_intensity_format, int_is_compressed, Int64(0), 0, :intensity)
                 spectra_metadata[j] = SpectrumMetadata(Int32(0), Int32(0), "", :sample, global_mode, mz_asset, int_asset)
             end
             break
@@ -431,8 +512,8 @@ function parse_imzml_spectrum_block(stream::IO, hIbd::IO, param_groups::Dict{Str
         
         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, 0.0, 0.0)
+            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, 0.0, 0.0)
+            int_asset = SpectrumAsset(default_intensity_format, int_is_compressed, Int64(0), 0, :intensity)
         else
             mz_info = mz_data[1]
             int_info = int_data[1]
@@ -446,9 +527,9 @@ function parse_imzml_spectrum_block(stream::IO, hIbd::IO, param_groups::Dict{Str
             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, 0.0, 0.0)
+                                    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, 0.0, 0.0)
+                                     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)
@@ -474,7 +555,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, use_mmap::Bool=true)
+function load_imzml_lazy(file_path::String; cache_size::Int=100)
     println("DEBUG: Checking for .imzML file at $file_path")
     if !isfile(file_path)
         throw(FileFormatError("Provided path is not a file: $(file_path)"))
@@ -488,118 +569,283 @@ function load_imzml_lazy(file_path::String; cache_size::Int=100, use_mmap::Bool=
 
     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
-        @debug "Parsing imzML header..."
+        # --- NEW: Parse all header information in a more efficient single pass ---
+        println("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
-        @debug "Image dimensions: $(width)x$(height), $num_spectra spectra."
+        println("DEBUG: Image dimensions: $(width)x$(height), $num_spectra spectra.")
 
-        # ... (format extraction logic stays the same) ...
+        # Extract default formats from the parsed param_groups
-        # [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; end
+            if group.Axis == 1
-        end
+                mz_group = group
-
+            elseif group.Axis == 2
-        default_mz_format = (mz_group !== nothing) ? mz_group.Format : Float64
+                int_group = group
-        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
 
-        # Check for compression status once
+        if mz_group === nothing || int_group === nothing
-        any_comp = mz_is_compressed || int_is_compressed
+            @warn "Could not find global definitions for m/z and intensity arrays. Using hardcoded defaults (Float64)."
-
+            default_mz_format = Float64
-        if !use_cache
+            default_intensity_format = Float64
-            @debug "Parsing spectrum block from XML (this may take time for large files)..."
+            mz_is_compressed = false
-            spectra_metadata = parse_imzml_spectrum_block(stream, ibd_handles[1], param_groups, width, height, num_spectra,
+            int_is_compressed = false
-                                                        default_mz_format, default_intensity_format, 
+            global_mode = UNKNOWN
-                                                        mz_is_compressed, int_is_compressed, global_mode)
+        else
-            
+            default_mz_format = mz_group.Format
-            @debug "Metadata parsing complete. Saving cache for next time..."
+            default_intensity_format = int_group.Format
-            # We create a temporary MSIData just for save_metadata_cache
+            mz_is_compressed = mz_group.Packed
-            tmp_source = ImzMLSource(ibd_handles, default_mz_format, default_intensity_format, mmap_data, any_comp)
+            int_is_compressed = int_group.Packed
-            tmp_msi = MSIData(tmp_source, spectra_metadata, instrument_meta, (width, height), nothing, cache_size)
+            global_mode = mz_group.Mode != UNKNOWN ? mz_group.Mode : int_group.Mode
-            save_metadata_cache(tmp_msi, cache_path)
         end
+
+        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,
+                                                          default_mz_format, default_intensity_format,
+                                                          mz_is_compressed, int_is_compressed, global_mode)
+
+        println("DEBUG: Metadata parsing complete.")
         
-        # Build coordinate map ...
+        # Build coordinate map for imzML files
+        println("DEBUG: Building 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.")
 
-        source = ImzMLSource(ibd_handles, default_mz_format, default_intensity_format, mmap_data, any_comp)
+        # --- NEW: Update acquisition mode based on spectrum parsing ---
-        @debug "Creating MSIData object."
+        acq_mode_symbol = if global_mode == CENTROID
-        msi_data = MSIData(source, spectra_metadata, instrument_meta, (width, height), coordinate_map, cache_size)
+            :centroid
+        elseif global_mode == PROFILE
+            :profile
+        else
+            :unknown
+        end
+        
+        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.
+        source = ImzMLSource(ts_hIbd, default_mz_format, default_intensity_format)
-        # The cache only stores binary offsets (SpectrumMetadataBinary). It does NOT populate
+        println("DEBUG: Creating MSIData object.")
-        # msi_data.spectrum_stats_df (TIC, BPI, BasePeakMZ, MinMZ, MaxMZ), which requires a
+        msi_data = MSIData(source, spectra_metadata, final_instrument_meta, (width, height), coordinate_map, cache_size)
-        # 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."
 
+        # Close the XML stream as it's no longer needed
         close(stream)
+
         return msi_data
 
     catch e
         close(stream)
-        # Check if handles exist before closing
+        close(ts_hIbd) # Ensure IBD handle is closed on error
-        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
+
 
 # =============================================================================
 #
@@ -623,7 +869,7 @@ This optimized version uses binary search for efficiency.
 # Returns
 - The intensity (`Float64`) of the peak if found, otherwise `0.0`.
 """
-@inline function find_mass(mz_array::AbstractVector{<:Real}, intensity_array::AbstractVector{<:Real}, 
+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)
@@ -684,91 +930,59 @@ 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
-        # Range check first (cheapest)
+        # NEW: Bloom filter check with discretization
-        # Access metadata directly if stats_df is missing to ensure zero-allocation
+        if bloom_filters !== nothing && !is_empty(bloom_filters[i])
-        @inbounds meta = data.spectra_metadata[i]
+            discretization_factor = 100.0
-        s_min, s_max = (stats_df !== nothing) ? (min_mzs[i], max_mzs[i]) : (meta.mz_asset.min_val, meta.mz_asset.max_val)
+            min_mass_int = round(Int, (mass - tolerance) * discretization_factor)
-        
+            max_mass_int = round(Int, (mass + tolerance) * discretization_factor)
-        if target_max < s_min || target_min > s_max
+            
-            continue
+            found = false
-        end
+            for mass_int in min_mass_int:max_mass_int
-
+                if mass_int in bloom_filters[i]
-        # Bloom filter rejection (very fast)
+                    found = true
-        if bloom_filters !== nothing
+                    break
-            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
-                @inbounds for mass_int in min_mass_int:max_mass_int
-                    if mass_int in bf
-                        found = true
-                        break
-                    end
                 end
-                !found && continue
+            end
+            
+            if !found
+                continue # Definitely not in this spectrum
             end
         end
 
-        candidate_count += 1
+        spec_min_mz = stats_df.MinMZ[i]
-        @inbounds candidate_indices[candidate_count] = i
+        spec_max_mz = stats_df.MaxMZ[i]
+        if target_max >= spec_min_mz && target_min <= spec_max_mz
+            push!(candidate_indices, i)
+        end
     end
 
-    # Use a view of the pre-allocated vector to avoid collect() allocations
+    println("Found $(length(candidate_indices)) candidate spectra (filtered from $(length(indices_to_check)) initial spectra)")
-    valid_candidates = view(candidate_indices, 1:candidate_count)
+
-    
     # 2. Iterate using the optimized, low-allocation iterator
-    # Use Atomic for thread-safe increment and avoid Ref-boxing
+    results_count = 0
-    results_count = Base.Threads.Atomic{Int}(0)
+    _iterate_spectra_fast(data, collect(candidate_indices)) do idx, mz_array, intensity_array
-    
+        meta = data.spectra_metadata[idx]
-    # Use let block to ensure closure captures are optimized (avoid boxing)
+        intensity = find_mass(mz_array, intensity_array, mass, tolerance)
-    let slice_matrix=slice_matrix, results_count=results_count, width=width, height=height, 
+        if intensity > 0.0
-        spectra_metadata=data.spectra_metadata, mass=mass, tolerance=tolerance
+            if 1 <= meta.x <= width && 1 <= meta.y <= height
-        
+                slice_matrix[meta.y, meta.x] = intensity
-        _iterate_spectra_fast(data, valid_candidates) do idx, mz_array, intensity_array
+                results_count += 1
-            @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
-                    @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
@@ -793,14 +1007,13 @@ 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 and allow binary search
+    # Sort masses to improve cache locality during search
     sorted_masses = sort(masses)
-    n_masses = length(sorted_masses)
 
-    # 1. Initialize a dictionary to hold the output slice matrices (using Float32 for 50% RAM savings)
+    # 1. Initialize a dictionary to hold the output slice matrices
-    slice_dict = Dict{Real, Matrix{Float32}}()
+    slice_dict = Dict{Real, Matrix{Float64}}()
     for mass in sorted_masses
-        slice_dict[mass] = zeros(Float32, height, width)
+        slice_dict[mass] = zeros(Float64, height, width)
     end
 
     local masked_indices::Union{Set{Int}, Nothing} = nothing
@@ -816,50 +1029,42 @@ function get_multiple_mz_slices(data::MSIData, masses::AbstractVector{<:Real}, t
         precompute_analytics(data)
     end
 
-    println("Filtering candidate spectra for $n_masses m/z values...")
+    println("Filtering candidate spectra for $(length(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. Optimized filtering: Iterate through spectra ONCE and check against all masses
+    # 3. Find all spectra that could contain *any* of the requested masses.
-    # This changes complexity from O(M*N) to O(N * log M) or O(N + M) depending on range overlap
+    for mass in sorted_masses
-    discretization_factor = 100.0
+        target_min = mass - tolerance
-    for i in indices_to_check
+        target_max = mass + tolerance
-        spec_min = stats_df.MinMZ[i]
+        for i in indices_to_check
-        spec_max = stats_df.MaxMZ[i]
+            # If already a candidate, no need to check again
-        
+            if i in candidate_indices
-        # Binary search to find masses that might overlap with this spectrum's range
+                continue
-        # target_min = mass - tolerance => mass = target_min + tolerance
+            end
-        # We need mass such that mass + tolerance >= spec_min  => mass >= spec_min - tolerance
+            # NEW: Bloom filter check with discretization
-        # 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
+                discretization_factor = 100.0
-                @inbounds for m_idx in m_start_idx:m_end_idx
+                min_mass_int = round(Int, (mass - tolerance) * discretization_factor)
-                    mass = sorted_masses[m_idx]
+                max_mass_int = round(Int, (mass + tolerance) * discretization_factor)
-                    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
-                    for mass_int in min_mass_int:max_mass_int
+                    if mass_int in bloom_filters[i]
-                        if mass_int in bloom_filters[i]
+                        found = true
-                            found_any = true
+                        break
-                            break
-                        end
                     end
-                    found_any && break
                 end
-                if found_any
+                
-                    push!(candidate_indices, i)
+                if !found
+                    continue # Definitely not in this spectrum
                 end
-            else
+            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
         end
@@ -868,40 +1073,29 @@ 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]
-        if isempty(mz_array)
+        # For this single spectrum, check all masses of interest
-            return
+        for mass in sorted_masses
-        end
+            # Check if this spectrum's range actually covers the current mass
-        
+            # This is a finer-grained check than the initial filtering
-        # Spectrum-level boundaries
+            if !isempty(mz_array) && (mass + tolerance) >= first(mz_array) && (mass - tolerance) <= last(mz_array)
-        spec_first = first(mz_array)
+                intensity = find_mass(mz_array, intensity_array, mass, tolerance)
-        spec_last = last(mz_array)
+                if intensity > 0.0
-        
+                    if 1 <= meta.x <= width && 1 <= meta.y <= height
-        # Find which of our target masses fall within this specific spectrum's actual range
+                        slice_dict[mass][meta.y, meta.x] = intensity
-        m_start_idx = searchsortedfirst(sorted_masses, spec_first - tolerance)
+                    end
-        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
-                    # 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 - replaces NaNs with 0.0 directly in Float32 matrices
+    # 5. Clean up and return
     for mass in sorted_masses
-        replace!(slice_dict[mass], NaN32 => 0.0f0)
+        replace!(slice_dict[mass], NaN => 0.0)
     end
     
-    println("Finished generating $n_masses slices in a single pass.")
+    println("Finished generating $(length(masses)) slices in a single pass.")
     return slice_dict
 end
 
@@ -1539,7 +1733,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{Real, Real}; 
+                               bounds::Tuple{Float64, Float64}; 
                                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, 0.0, 0.0)
+    return SpectrumAsset(data_format, compression_flag, binary_offset, encoded_length, axis)
 end
 
 # This function is updated to return the generic SpectrumMetadata struct
@@ -394,48 +394,38 @@ 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(primary_handle)
+        instrument_meta = parse_instrument_metadata_mzml(ts_stream.handle)
 
-        seekstart(primary_handle) # Reset stream after header 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("------------------------------------")
+
+        seekstart(ts_stream.handle) # Reset stream after header parsing
 
         println("DEBUG: Finding index offset...")
-        index_offset = find_index_offset(primary_handle)
+        index_offset = find_index_offset(ts_stream.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(primary_handle, index_offset)
+        seek(ts_stream.handle, index_offset)
         
         println("DEBUG: Searching for '<index name=\"spectrum\">'.")
-        if find_tag(primary_handle, r"<index\s+name=\"spectrum\"") === nothing
+        if find_tag(ts_stream.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(primary_handle)
+        spectrum_offsets = parse_offset_list(ts_stream.handle)
         if isempty(spectrum_offsets)
             throw(FileFormatError("No spectrum offsets found."))
         end
@@ -443,25 +433,31 @@ 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(primary_handle, spectrum_offsets[i])
+            spectra_metadata[i] = parse_spectrum_metadata(ts_stream.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.")
 
-        # Inferred global formats from first spectrum
+        # Assuming uniform data formats, take from the 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")
         
-        # Determine overall acquisition mode ...
+        # --- NEW: Determine overall acquisition mode ---
-        num_centroid = count(m -> m.mode == CENTROID, spectra_metadata)
+        modes = [meta.mode for meta in spectra_metadata]
-        num_profile = count(m -> m.mode == PROFILE, spectra_metadata)
+        num_centroid = count(m -> m == CENTROID, modes)
+        num_profile = count(m -> m == PROFILE, modes)
         
         acq_mode_symbol = if num_centroid > 0 && num_profile == 0
             :centroid
@@ -472,28 +468,26 @@ 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, 
+            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
+            instrument_meta.vendor_preprocessing_steps # Add this new field
         )
 
-        source = MzMLSource(mzml_handles, mz_format, intensity_format, mmap_data)
+        source = MzMLSource(ts_stream, mz_format, intensity_format)
         println("DEBUG: Creating MSIData object.")
         return MSIData(source, spectra_metadata, final_instrument_meta, (0, 0), nothing, cache_size)
 
     catch e
-        # Close all handles in the pool if initialization fails
+        close(ts_stream) # Ensure stream is closed on error
-        for h in mzml_handles
-            isopen(h) && close(h)
-        end
         rethrow(e)
     end
 end

start_MSI_GUI.jl

@@ -19,52 +19,6 @@ 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,16 +49,12 @@ const BENCHMARK_CASES = [
     # BenchmarkCase("/path/to/your/small_file.imzML", 309.06, 0.1),
     # BenchmarkCase("/path/to/your/medium_file.imzML", 896.0, 1.0),
     # BenchmarkCase("/path/to/your/large_file.imzML", 100.0, 0.1),
-    # BenchmarkCase("/home/pixel/Documents/Cinvestav_2025/Analisis/imzML_LA-ESI/180817_NEG_Thaliana_Leaf_bottom_1_0841.imzML",116.07,0.1, "Thaliana Leaf"),
+    BenchmarkCase("/home/pixel/Documents/Cinvestav_2025/Analisis/imzML_LA-ESI/180817_NEG_Thaliana_Leaf_bottom_1_0841.imzML",116.07,0.1, "Thaliana Leaf"),
-    BenchmarkCase("/home/pixel/Documents/Cinvestav_2025/Analisis/Chilli/ltpmsi-chilli.imzML",420,0.1, "Chilli Pepper"), # Chilli
+    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/40TopL,10TopR,30BottomL,20BottomR/40TopL,10TopR,30BottomL,20BottomR-centroid.imzML",885.5,0.1, "Colon Cancer Human"), # Human Cancer
+    BenchmarkCase("/home/pixel/Documents/Cinvestav_2025/Analisis/imzML_DESI/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_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 2GB"), # Mouse Stomach
+    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_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

@@ -1,215 +0,0 @@
-# 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

@@ -1,48 +0,0 @@
-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,9 +16,8 @@ using MSI_src
 const TEST_MZML_FILE = ""
 # const TEST_IMZML_FILE = "/home/pixel/Documents/Cinvestav_2025/Analisis/CE4_BF_R1/CE4_BF_R1.imzML"
 #const TEST_IMZML_FILE = "/home/pixel/Documents/Cinvestav_2025/Analisis/salida/Stomach_DHB_uncompressed.imzML"
-# const TEST_IMZML_FILE = "/home/pixel/Documents/Cinvestav_2025/Analisis/Thricoderma_etc/Imaging_interaccion_trichoderma_vs_streptomyces.imzML"
+const TEST_IMZML_FILE = "/home/pixel/Documents/Cinvestav_2025/Analisis/Thricoderma_etc/Imaging_interaccion_trichoderma_vs_streptomyces.imzML"
 #const MASK_ROUTE = "/home/pixel/Documents/Cinvestav_2025/JuliaMSI/public/css/masks/Stomach_DHB_uncompressed.png"
-const TEST_IMZML_FILE = "/home/pixel/Documents/Cinvestav_2025/Analisis/imzML_AP_SMALDI/HR2MSImouseurinarybladderS096.imzML"
 const MASK_ROUTE = ""
 
 #=

test/run_preprocessing.jl

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

test/run_tests.jl

@@ -31,10 +31,9 @@ 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 ---
@@ -59,17 +58,16 @@ 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.2
+const MZ_TOLERANCE = 0.1
 
 # Coordinates to plot a specific spectrum from imzML
 const COORDS_TO_PLOT = (50, 50) # Example coordinates (X, Y)

test/runtests.jl

@@ -1,58 +0,0 @@
-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

@@ -1,124 +0,0 @@
-#!/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()