cleaned the repository from dead code, simplified some functions, users can now create a system image of the JuliaMSI repo that runs the core dependencies faster without compiler time, alongside a docker environment for running scripts, the instructions are on the readme

.github/workflows/ci.yml

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

Dockerfile.headless

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

Project.toml

@@ -1,4 +1,5 @@
 name = "MSI_src"
+uuid = "8f395f85-3b9c-4f7f-bf7e-12345678abcd"
 authors = ["JJSA"]
 version = "0.1.0"
 
@@ -47,6 +48,7 @@ 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]
@@ -92,5 +94,12 @@ Setfield = "1.1"
 Statistics = "1.11"
 StatsBase = "0.34"
 StipplePlotly = "0.13"
+Test = "1.11.0"
 UUIDs = "1.11"
 julia = "1.11, 1.12"
+
+[extras]
+Test = "8dfed614-e60c-5e00-aba5-a33c2d43236e"
+
+[targets]
+test = ["Test"]

README.md

@@ -35,35 +35,44 @@ 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 the .so file by running the next build script in your directory:
+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–10 minutes as it merges all dependencies and JuliaMSI functions into a single binary. The resulting .so/.dll file will be around **300MB–600MB** because it contains the pre-compiled machine code for your entire environment.
+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.
-A sysimage built on Linux (.so) will not work on Windows.
+*Note: A sysimage built on Linux (.so) will not work on Windows. You must run the build script once on each target operating system.*
-You must run the build_sysimage.jl script once on each target operating system.
 
-Once MSI_sysimage.so is created in your directory, adapt your command like this:
+Once `sys_msi_gui.so` is created in your directory, adapt your launch command:
 ```bash
 julia --project=. -e 'using Pkg; Pkg.precompile()'
-# This will find MSI_sysimage.so, .dll, or .dylib automatically:
+# Adjust the .so extension to .dll if on Windows or .dylib if on macOS
-julia --threads auto --project=. --sysimage MSI_sysimage* start_MSI_GUI.jl
+julia --threads auto --project=. --sysimage sys_msi_gui.so start_MSI_GUI.jl
 ```
-The command above loads a pre-compiled version of your environment (note, you have to use either the sysimage command or the normal start command), speeding up the booting delay. Your GUI or scripts using JuliaMSI should start almost instantly and process files significantly faster than before.
-Always run the build script **in the same project root** where you intend to run the app.
 
-In scripts like test/run_tests.jl, you no longer need to manually include("../src/MSI_src.jl"). You can simply treat it like a globally installed package to load the precompiled binary:
+### 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:
-```julia
-using MSI_src
-# ... rest of your script
-```
-### Run with the --sysimage flag
-Launch your tests using the same flag to bypass all compilation overhead:
 
 ```bash
-julia --threads auto --project=. --sysimage MSI_sysimage.so test/run_tests.jl
+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

app.jl.html

@@ -546,13 +546,19 @@
           </q-list>
 
           <!-- Pipeline Controls -->
-          <div class="row justify-end items-center q-mt-md">
+          <div class="row q-mt-md justify-between items-center text-caption text-grey">
-            <q-btn class="q-ma-sm" icon="get_app" v-on:click="export_params_btn=true" label="Export Params" outline
+            <div>
-              :disable="is_processing" />
+               <q-icon name="info" class="q-mr-xs" /> Feature matrices are exported in highly optimized MatrixMarket Coordinate format (<b>.mtx</b>).
-            <q-btn class="q-ma-sm" icon="upload_file" v-on:click="import_params_btn=true" label="Import Params" outline
+               <br/><i>Python:</i> <code>scipy.io.mmread('filtered_data.mtx')</code> | <i>R:</i> <code>Matrix::readMM('filtered_data.mtx')</code>
-              :disable="is_processing" />
+            </div>
-            <q-btn :loading="is_processing" class="q-ma-sm btn-style" icon="play_arrow"
+            <div>
-              v-on:click="run_full_pipeline=true" padding="lg" label="Run Pipeline" :disable="is_processing" />
+              <q-btn class="q-ma-sm" icon="get_app" v-on:click="export_params_btn=true" label="Export Params" outline
+                :disable="is_processing" />
+              <q-btn class="q-ma-sm" icon="upload_file" v-on:click="import_params_btn=true" label="Import Params" outline
+                :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,27 +1,35 @@
 # build_sysimage.jl
 using PackageCompiler
 
+headless = "--headless" in ARGS
+
 println("--- Starting Custom System Image Build ---")
-println("This process can take 5-10 minutes.")
+println("Mode: ", headless ? "Headless (Scripting/Docker)" : "Full GUI (Genie App)")
+println("This process can take 5-15 minutes.")
 
-# Define the image path with OS-specific extension
 ext = Sys.iswindows() ? ".dll" : Sys.isapple() ? ".dylib" : ".so"
-sysimage_path = "MSI_sysimage" * ext
+sysimage_name = headless ? "sys_msi_headless" * ext : "sys_msi_gui" * ext
+sysimage_path = joinpath(@__DIR__, sysimage_name)
 
-# List dependencies to include in the sysimage for maximum stability
+packages_to_include = Symbol[
-# We have removed graphical heavy-lifters (CairoMakie, PlotlyBase) to ensure 
+    :MSI_src, # <--- Bake our core library
-# the build completes successfully on systems with standard RAM.
-packages_to_include = [
-    :Genie,
-    :GenieFramework,
-    :PlotlyBase,
-    :Libz,
     :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,
@@ -31,5 +39,11 @@ create_sysimage(
 )
 
 println("--- System Image Build Successful! ---")
-println("To start Julia with this image, use:")
+println("Created: $(sysimage_path)")
-println("julia --threads auto --project=. --sysimage $(sysimage_path) start_MSI_GUI.jl")
+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

precompile_script.jl

@@ -133,6 +133,22 @@ try
                 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))
+                    ],
+                    adaptive_binning=true,
+                    bin_tolerance_ppm=50.0
+                )
+                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)

src/BloomFilters.jl

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

src/Common.jl

@@ -31,76 +31,7 @@ function posix_madvise(buffer::AbstractArray, advice::Integer)
 end
 
 # --- Buffer Pooling ---
-"""
-    BufferPool
 
-A thread-safe pool of `Vector{UInt8}` buffers, categorized by their size.
-This helps reduce memory allocations and garbage collection overhead by reusing buffers.
-
-# Arguments:
-- `lock`: A `ReentrantLock` to ensure thread-safe access to the pool.
-- `buffers`: A dictionary mapping buffer size (in bytes) to a list of available buffers of that size.
-"""
-mutable struct BufferPool
-    lock::ReentrantLock
-    buffers::Dict{Int, Vector{Vector{UInt8}}}
-
-    function BufferPool()
-        new(ReentrantLock(), Dict{Int, Vector{Vector{UInt8}}}())
-    end
-end
-
-"""
-    get_buffer!(pool::BufferPool, size::Int) -> Vector{UInt8}
-
-Retrieves a `Vector{UInt8}` buffer of at least `size` bytes from the pool.
-If no suitable buffer is available, a new one is allocated.
-
-# Arguments:
-- `pool`: The `BufferPool` to retrieve the buffer from.
-- `size`: The minimum desired size of the buffer in bytes.
-
-# Returns:
-- A `Vector{UInt8}` buffer.
-"""
-function get_buffer!(pool::BufferPool, size::Int)::Vector{UInt8}
-    lock(pool.lock) do
-        # Find an existing buffer that is large enough
-        for (buffer_size, buffer_list) in pool.buffers
-            if buffer_size >= size && !isempty(buffer_list)
-                buffer = pop!(buffer_list)
-                # Resize if necessary (and if it's significantly larger than needed)
-                if length(buffer) > size * 2 # Heuristic: if buffer is more than twice the requested size
-                    return Vector{UInt8}(undef, size) # Allocate new smaller buffer
-                else
-                    return buffer
-                end
-            end
-        end
-        # No suitable buffer found, allocate a new one
-        return Vector{UInt8}(undef, size)
-    end
-end
-
-"""
-    release_buffer!(pool::BufferPool, buffer::Vector{UInt8})
-
-Returns a `Vector{UInt8}` buffer to the pool for future reuse.
-
-# Arguments:
-- `pool`: The `BufferPool` to return the buffer to.
-- `buffer`: The `Vector{UInt8}` buffer to release.
-"""
-function release_buffer!(pool::BufferPool, buffer::Vector{UInt8})
-    lock(pool.lock) do
-        buffer_size = length(buffer)
-        if !haskey(pool.buffers, buffer_size)
-            pool.buffers[buffer_size] = Vector{Vector{UInt8}}()
-        end
-        push!(pool.buffers[buffer_size], buffer)
-    end
-    return nothing
-end
 
 
 """
@@ -140,11 +71,16 @@ If no suitable buffer is available, a new one is allocated.
 - A `Vector{UInt8}` buffer.
 """
 function get_buffer!(pool::SimpleBufferPool, size::Int)::Vector{UInt8}
-    lock(pool.lock) do
+    buf = lock(pool.lock) do
         # Check for existing buffers of exact size first
         if haskey(pool.buffers, size) && !isempty(pool.buffers[size])
             return pop!(pool.buffers[size])
         end
+        return nothing
+    end
+    
+    if buf !== nothing
+        return buf
     end
     
     # No suitable buffer found, allocate new one (outside lock to reduce contention)

src/FusedPipeline.jl

@@ -1,93 +0,0 @@
-# src/FusedPipeline.jl
-# This file defines a high-performance in-place preprocessing pipeline
-# that minimizes allocations by using reused buffers from ResourcePool.
-
-using .MSI_src # Ensure it can see the module's contents if needed
-
-"""
-    SpectralPipeline
-
-Holds a sequence of preprocessing steps and the necessary buffers to execute them in-place.
-"""
-struct SpectralPipeline
-    steps::Vector{AbstractPreprocessingStep}
-end
-
-"""
-    apply_pipeline!(mz::Vector{Float64}, intensity::Vector{Float64}, pipeline::SpectralPipeline; data::MSIData)
-
-Applies all steps in the pipeline to the spectrum arrays in-place.
-Uses internal buffers from data.resource_pool where needed.
-"""
-function apply_pipeline!(mz::Vector{Float64}, intensity::Vector{Float64}, pipeline::SpectralPipeline, data::MSIData)
-    # Process each step in sequence
-    for step in pipeline.steps
-        apply_step!(mz, intensity, step, data)
-    end
-    return mz, intensity
-end
-
-# --- Basic in-place implementations of core preprocessing steps ---
-
-function apply_step!(mz, int, step::Normalization, data)
-    if step.method === :tic
-        s = sum(int)
-        if s > 0
-            int ./= s
-        end
-    elseif step.method === :median
-        m = median(int)
-        if m > 0
-            int ./= m
-        end
-    end
-    return int
-end
-
-function apply_step!(mz, int, step::Smoothing, data)
-    if step.method === :savitzky_golay
-        # SavitzkyGolay.savitzky_golay currently allocates, but we can't easily fix that here
-        # without refactoring the library. However, we can use a pooled vector for its output
-        # then copy back to intensity.
-        # [Wait: For now we'll call smoothed_y = smooth_spectrum_core(int, ...)]
-        # We'll use a resource from the pool to avoid fresh allocation
-        temp_buf = acquire(data.resource_pool)
-        resize!(temp_buf, length(int))
-        
-        # Call existing core which returns a new vector, unfortunately
-        # But we'll copy it back to 'int' to maintain in-place pipeline
-        smoothed = smooth_spectrum_core(int; method=step.method, window=step.window, order=step.order)
-        copyto!(int, smoothed)
-        
-        release!(data.resource_pool, temp_buf)
-    end
-    return int
-end
-
-function apply_step!(mz, int, step::BaselineCorrection, data)
-    if step.method === :snip
-        # SNIP is easy to make in-place!
-        iterations = (step.iterations === nothing) ? 100 : step.iterations
-        snip_baseline_inplace!(int, iterations)
-    end
-    return int
-end
-
-"""
-    snip_baseline_inplace!(y, iterations)
-
-In-place implementation of the Sensitive Nonlinear Iterative Peak clipping algorithm.
-"""
-function snip_baseline_inplace!(y::Vector{Float64}, iterations::Int)
-    n = length(y)
-    n < 3 && return y
-    
-    # We still need one temporary buffer for the SNIP iteration to read from the previous state
-    # Actually, we can just return the baseline and subtract it, but to BE in-place,
-    # we need a temporary to hold the baseline during calculation.
-    
-    # For now, we'll use the existing _snip_baseline_impl and subtract
-    baseline = _snip_baseline_impl(y, iterations=iterations)
-    y .-= baseline
-    return y
-end

src/MSIData.jl

@@ -10,100 +10,6 @@ using Base64, Libz, Serialization, Printf, DataFrames, Base.Threads, StatsBase,
 
 const FILE_HANDLE_LOCK = ReentrantLock()
 
-"""
-    ThreadSafeFileHandle
-
-Wraps a file handle with thread-safe operations.
-# Fields
-
-- `path::String`: The path to the file.
-- `handle::IO`: The file handle.
-- `lock::ReentrantLock`: The lock for thread-safe operations.
-"""
-mutable struct ThreadSafeFileHandle
-    path::String
-    handle::IO
-    lock::ReentrantLock
-end
-
-"""
-    ThreadSafeFileHandle(path::String, mode::String="r")
-
-Wraps a file handle with thread-safe operations.
-
-# Arguments
-
-- `path::String`: The path to the file.
-- `mode::String`: The mode to open the file in.
-
-# Returns
-
-- A `ThreadSafeFileHandle`.
-"""
-function ThreadSafeFileHandle(path::String, mode::String="r")
-    handle = lock(FILE_HANDLE_LOCK) do
-        open(path, mode)
-    end
-    return ThreadSafeFileHandle(path, handle, ReentrantLock())
-end
-
-"""
-    read_at!(tsfh::ThreadSafeFileHandle, a::AbstractArray, pos::Integer)
-
-Atomically seeks to a position and reads data into an array. This is the thread-safe
-way to read from a specific offset in the file.
-"""
-function read_at!(tsfh::ThreadSafeFileHandle, a::AbstractArray, pos::Integer)
-    if pos < 0
-        throw(ArgumentError("Invalid seek position: $pos"))
-    end
-    lock(tsfh.lock) do
-        # Also check against file size if possible, though filesize() might be expensive to call repeatedly
-        # Let's rely on seek throwing if it goes way out of bounds, but catch the negative case which is definitely an error.
-        seek(tsfh.handle, pos)
-        read!(tsfh.handle, a)
-    end
-end
-
-function Base.read(tsfh::ThreadSafeFileHandle, n::Integer)
-    lock(tsfh.lock) do
-        return read(tsfh.handle, n)
-    end
-end
-
-function Base.close(tsfh::ThreadSafeFileHandle)
-    close(tsfh.handle)
-end
-
-function Base.isopen(tsfh::ThreadSafeFileHandle)
-    isopen(tsfh.handle)
-end
-
-function Base.position(tsfh::ThreadSafeFileHandle)
-    lock(tsfh.lock) do
-        return position(tsfh.handle)
-    end
-end
-
-function Base.filesize(tsfh::ThreadSafeFileHandle)
-    lock(tsfh.lock) do
-        return filesize(tsfh.handle)
-    end
-end
-
-# Overload Base.seek for ThreadSafeFileHandle
-function Base.seek(tsfh::ThreadSafeFileHandle, pos::Integer)
-    lock(tsfh.lock) do
-        seek(tsfh.handle, pos)
-    end
-end
-
-# Overload Base.read! for ThreadSafeFileHandle
-function Base.read!(tsfh::ThreadSafeFileHandle, a::AbstractArray)
-    lock(tsfh.lock) do
-        read!(tsfh.handle, a)
-    end
-end
 
 # Abstract type for different data sources (e.g., mzML, imzML)
 # This allows dispatching to the correct binary reading logic.
@@ -331,7 +237,7 @@ mutable struct MSIData
     global_min_mz::Base.Threads.Atomic{Float64}
     global_max_mz::Base.Threads.Atomic{Float64}
     spectrum_stats_df::Union{DataFrame, Nothing}
-    bloom_filters::Union{Vector{<:BloomFilter}, Nothing}
+    bloom_filters::Union{Vector{Union{BloomFilter{Int}, Nothing}}, Nothing}
     analytics_ready::AtomicFlag
     preprocessing_hints::Union{Dict{Symbol, Any}, Nothing} # For auto-determined parameters
 
@@ -487,65 +393,7 @@ function Base.close(data::MSIData)
     println("MSIData object closed and resources released.")
 end
 
-"""
-    warm_cache(data::MSIData, indices::AbstractVector{Int})
 
-Pre-loads a specified list of spectra into the cache. This is useful for warming up
-the cache with frequently accessed spectra to improve performance for subsequent
-interactive analysis.
-
-# Arguments
-
-- `data::MSIData`: The MSIData object.
-- `indices::AbstractVector{Int}`: The indices of the spectra to warm the cache with.
-
-# Returns
-- `nothing`
-"""
-function warm_cache(data::MSIData, indices::AbstractVector{Int})
-    println("Warming cache with $(length(indices)) spectra...")
-    for idx in indices
-        # Calling GetSpectrum will load the data and place it in the cache
-        GetSpectrum(data, idx)
-    end
-    println("Cache warming complete.")
-end
-
-"""
-    get_memory_usage(data::MSIData)
-
-Calculates and returns a summary of the memory currently used by the MSIData object,
-including the spectrum cache and metadata.
-
-# Arguments
-
-- `data::MSIData`: The MSIData object.
-
-# Returns
-- `(total_mb, cache_mb, metadata_mb, coordinate_map_mb)`: A named tuple with memory usage in MB for different components.
-"""
-function get_memory_usage(data::MSIData)
-    cache_bytes = 0
-    lock(data.cache_lock) do
-        # This is an approximation; `sizeof` on a dictionary is not fully representative,
-        # but it's much better to sum the size of the actual vectors.
-        for (mz, intensity) in values(data.cache)
-            cache_bytes += sizeof(mz) + sizeof(intensity)
-        end
-    end
-
-    metadata_bytes = sizeof(data.spectra_metadata)
-    coordinate_map_bytes = data.coordinate_map === nothing ? 0 : sizeof(data.coordinate_map)
-
-    total_bytes = cache_bytes + metadata_bytes + coordinate_map_bytes
-
-    return (
-        total_mb = total_bytes / 1024^2,
-        cache_mb = cache_bytes / 1024^2,
-        metadata_mb = metadata_bytes / 1024^2,
-        coordinate_map_mb = coordinate_map_bytes / 1024^2
-    )
-end
 
 """
     get_masked_spectrum_indices(data::MSIData, mask_matrix::BitMatrix) -> Set{Int}
@@ -648,26 +496,7 @@ function read_binary_vector(data::MSIData, io::IO, asset::SpectrumAsset)
     return out_array
 end
 
-"""
-    read_binary_vector(data::MSIData, ts_handle::ThreadSafeFileHandle, asset::SpectrumAsset)
 
-Reads a single spectrum's m/z and intensity arrays directly from the `.ibd`
-binary file for an `.imzML` dataset. It handles reading the raw binary data
-and converting it from little-endian to the host's native byte order.
-
-# Arguments
-- `data`: The `MSIData` object.
-- `ts_handle`: The `ThreadSafeFileHandle` containing the file handle and data formats.
-- `asset`: The `SpectrumAsset` containing metadata for the array.
-
-# Returns
-- A `Vector` of the appropriate type containing the decoded data.
-"""
-function read_binary_vector(data::MSIData, ts_handle::ThreadSafeFileHandle, asset::SpectrumAsset)
-    lock(ts_handle.lock) do
-        return read_binary_vector(data, ts_handle.handle, asset)
-    end
-end
 
 # Overload for different source types
 """
@@ -747,9 +576,10 @@ This is a high-level function that orchestrates the reading process by calling
 - A tuple `(mz, intensity)` containing the two requested data arrays.
 """
 function read_spectrum_from_disk(data::MSIData, source::MzMLSource, meta::SpectrumMetadata)
+    handle = get_handle(source)
     # For mzML, data is Base64 encoded within the XML
-    mz = read_binary_vector(data, source.file_handle, meta.mz_asset)
+    mz = read_binary_vector(data, handle, meta.mz_asset)
-    intensity = read_binary_vector(data, source.file_handle, meta.int_asset)
+    intensity = read_binary_vector(data, handle, meta.int_asset)
     
     mz_f64, intensity_f64 = Float64.(mz), Float64.(intensity)
 
@@ -1064,8 +894,9 @@ function precompute_analytics(msi_data::MSIData)
         max_mzs = Vector{Float64}(undef, num_spectra)
         modes = Vector{SpectrumMode}(undef, num_spectra)
         
-        # Don't precompute Bloom filters by default - create them lazily
+        # Precompute Bloom filters during this pass
-        bloom_filters = nothing
+        bloom_filters = Vector{Union{BloomFilter{Int}, Nothing}}(undef, num_spectra)
+        fill!(bloom_filters, nothing)
 
         # Process in chunks to reduce memory pressure
         chunk_size = min(10000, num_spectra)
@@ -1105,6 +936,7 @@ function precompute_analytics(msi_data::MSIData)
                 bpis[idx] = max_int
                 bp_mzs[idx] = mz[max_idx]
                 num_points[idx] = length(mz)
+                bloom_filters[idx] = create_bloom_filter_for_spectrum(mz)
             end
             
             # Force GC between chunks to control memory
@@ -1163,17 +995,17 @@ end
 
 
 
+
 """
-    get_total_spectrum_imzml(msi_data::MSIData; num_bins::Int=2000)
+    get_total_spectrum_core(msi_data::MSIData; num_bins::Int=2000, masked_indices::Union{Set{Int}, Nothing}=nothing)
 
-Internal function to calculate the total spectrum for an `.imzML` file.
+Internal function to calculate the total spectrum.
 
-It uses a fast, two-pass approach:
+It uses a two-pass approach:
-1. The first pass finds the global m/z range across all spectra.
+1. The first pass finds the global m/z range by iterating through all spectra.
 2. The second pass sums intensities into a pre-defined number of bins.
 
-This function is highly optimized for `.imzML` by leveraging direct binary
+This function is called by `get_total_spectrum`.
-reading and optimized binning logic. It is called by `get_total_spectrum`.
 
 # Arguments
 
@@ -1182,16 +1014,10 @@ reading and optimized binning logic. It is called by `get_total_spectrum`.
 - `masked_indices::Union{Set{Int}, Nothing}`: The indices of the spectra to mask.
 
 # Returns
-- `(mz_range, total_spectrum)`: A tuple containing the m/z range and the total spectrum.
+- `(mz_range, total_spectrum, num_spectra_processed)`: A tuple containing the m/z range, total spectrum, and spectrum count.
-
-# Example
-
-```julia
-get_total_spectrum_imzml(msi_data)
-```
 """
-function get_total_spectrum_imzml(msi_data::MSIData; num_bins::Int=2000, masked_indices::Union{Set{Int}, Nothing}=nothing)
+function get_total_spectrum_core(msi_data::MSIData; num_bins::Int=2000, masked_indices::Union{Set{Int}, Nothing}=nothing)
-    println("Calculating total spectrum for imzML (2-pass method)...")
+    println("Calculating total spectrum (2-pass method)...")
     total_start_time = time_ns()
 
     local global_min_mz, global_max_mz
@@ -1201,7 +1027,7 @@ function get_total_spectrum_imzml(msi_data::MSIData; num_bins::Int=2000, masked_
         global_min_mz = Base.Threads.atomic_add!(msi_data.global_min_mz, 0.0)
         global_max_mz = Base.Threads.atomic_add!(msi_data.global_max_mz, 0.0)
     else
-        # 1. First Pass: Find the global m/z range by reading the fast .ibd file
+        # 1. First Pass: Find the global m/z range
         pass1_start_time = time_ns()
         println("  Pass 1: Finding global m/z range...")
         g_min_mz = Inf
@@ -1214,8 +1040,9 @@ function get_total_spectrum_imzml(msi_data::MSIData; num_bins::Int=2000, masked_
             end
         end
         pass1_duration = (time_ns() - pass1_start_time) / 1e9
+        
         if !isfinite(g_min_mz)
-            @warn "Could not determine a valid m/z range for imzML. All spectra might be empty."
+            @warn "Could not determine a valid m/z range. All spectra might be empty."
             return (Float64[], Float64[], 0)
         end
         
@@ -1226,6 +1053,7 @@ function get_total_spectrum_imzml(msi_data::MSIData; num_bins::Int=2000, masked_
         println("  Global m/z range found and cached: [$(global_min_mz), $(global_max_mz)]")
         @printf "  (Pass 1 took %.2f seconds)\n" pass1_duration
     end
+    
     # 2. Define Bins and precompute constants
     mz_bins = range(global_min_mz, stop=global_max_mz, length=num_bins)
     bin_step = step(mz_bins)
@@ -1272,127 +1100,12 @@ function get_total_spectrum_imzml(msi_data::MSIData; num_bins::Int=2000, masked_
     pass2_duration = (time_ns() - pass2_start_time) / 1e9
 
     total_duration = (time_ns() - total_start_time) / 1e9
-    println("\n--- imzML Profiling (Post-Optimization) ---")
+    println("\n--- Profiling (Post-Optimization) ---")
     @printf "  Pass 2 (I/O + Binning):   %.2f seconds\n" pass2_duration
     @printf "  Total Function Time:      %.2f seconds\n" total_duration
     println("----------------------------------------\n")
 
-    println("Total spectrum calculation complete for imzML.")
+    println("Total spectrum calculation complete.")
-    return (collect(mz_bins), intensity_sum, num_spectra_processed)
-end
-
-"""
-    get_total_spectrum_mzml(msi_data::MSIData; num_bins::Int=2000)
-
-Internal function to calculate the total spectrum for an `.mzML` file.
-
-It uses a two-pass approach analogous to the `imzML` implementation:
-1. The first pass finds the global m/z range by iterating through all spectra.
-2. The second pass sums intensities into a pre-defined number of bins.
-
-This function is called by `get_total_spectrum`.
-
-# Arguments
-
-- `msi_data::MSIData`: The MSIData object.
-- `num_bins::Int`: The number of bins to use for the total spectrum.
-- `masked_indices::Union{Set{Int}, Nothing}`: The indices of the spectra to mask.
-
-# Returns
-- `(mz_range, total_spectrum)`: A tuple containing the m/z range and the total spectrum.
-
-# Example
-
-```julia
-get_total_spectrum_mzml(msi_data)
-```
-"""
-function get_total_spectrum_mzml(msi_data::MSIData; num_bins::Int=2000, masked_indices::Union{Set{Int}, Nothing}=nothing)
-    println("Calculating total spectrum for mzML (2-pass method)...")
-    total_start_time = time_ns()
-
-    local global_min_mz, global_max_mz
-
-    if Base.Threads.atomic_add!(msi_data.global_min_mz, 0.0) !== Inf
-        println("  Using pre-computed m/z range.")
-        global_min_mz = Base.Threads.atomic_add!(msi_data.global_min_mz, 0.0)
-        global_max_mz = Base.Threads.atomic_add!(msi_data.global_max_mz, 0.0)
-    else
-        # --- Pass 1: Find m/z range and cache it ---
-        pass1_start_time = time_ns()
-        println("  Pass 1: Finding global m/z range...")
-        g_min_mz = Inf
-        g_max_mz = -Inf
-        _iterate_spectra_fast(msi_data, masked_indices === nothing ? nothing : collect(masked_indices)) do idx, mz, intensity
-            if !isempty(mz)
-                local_min, local_max = extrema(mz)
-                g_min_mz = min(g_min_mz, local_min)
-                g_max_mz = max(g_max_mz, local_max)
-            end
-        end
-        pass1_duration = (time_ns() - pass1_start_time) / 1e9
-
-        if !isfinite(g_min_mz)
-            @warn "Could not determine a valid m/z range for mzML. All spectra might be empty."
-            return (Float64[], Float64[], 0)
-        end
-
-        # Use and cache the result
-        global_min_mz = g_min_mz
-        global_max_mz = g_max_mz
-        set_global_mz_range!(msi_data, global_min_mz, global_max_mz)
-
-        println("  Global m/z range found and cached: [$(global_min_mz), $(global_max_mz)]")
-        @printf "  (Pass 1 took %.2f seconds)\n" pass1_duration
-    end
-
-    # --- Pass 2: Bin intensities ---
-    pass2_start_time = time_ns()
-    println("  Pass 2: Summing intensities into $num_bins bins...")
-    
-    mz_bins = range(global_min_mz, stop=global_max_mz, length=num_bins)
-    bin_step = step(mz_bins)
-    inv_bin_step = 1.0 / bin_step  # Precompute reciprocal
-    min_mz = global_min_mz
-
-    # Initialize thread-local intensity sums
-    thread_local_intensity_sums = [zeros(Float64, num_bins) for _ in 1:Base.Threads.nthreads()]
-    thread_local_spectra_processed = [0 for _ in 1:Base.Threads.nthreads()]
-
-    _iterate_spectra_fast(msi_data, masked_indices === nothing ? nothing : collect(masked_indices)) do idx, mz, intensity
-        thread_id = Base.Threads.threadid()
-        thread_local_spectra_processed[thread_id] += 1
-        if isempty(mz)
-            return
-        end
-        
-        @inbounds @simd for i in eachindex(mz)
-            # Calculate raw bin index
-            raw_index = (mz[i] - min_mz) * inv_bin_step + 1.0
-            bin_index = trunc(Int, raw_index)
-            
-            # Branchless version using clamp
-            final_index = clamp(bin_index, 1, num_bins)
-            thread_local_intensity_sums[thread_id][final_index] += intensity[i]
-        end
-    end
-
-    # Combine thread-local results
-    intensity_sum = zeros(Float64, num_bins)
-    num_spectra_processed = 0
-    for i in 1:Base.Threads.nthreads()
-        intensity_sum .+= thread_local_intensity_sums[i]
-        num_spectra_processed += thread_local_spectra_processed[i]
-    end
-    pass2_duration = (time_ns() - pass2_start_time) / 1e9
-
-    total_duration = (time_ns() - total_start_time) / 1e9
-    println("\n--- mzML Profiling (Post-Optimization) ---")
-    @printf "  Pass 2 (I/O + Binning):   %.2f seconds\n" pass2_duration
-    @printf "  Total Function Time:      %.2f seconds\n" total_duration
-    println("----------------------------------------\n")
-
-    println("Total spectrum calculation complete for mzML.")
     return (collect(mz_bins), intensity_sum, num_spectra_processed)
 end
 
@@ -1428,11 +1141,7 @@ function get_total_spectrum(msi_data::MSIData; num_bins::Int=2000, mask_path::Un
         println("Calculating total spectrum for masked region (mask from: $(mask_path))")
     end
 
-    if msi_data.source isa ImzMLSource
+    return get_total_spectrum_core(msi_data, num_bins=num_bins, masked_indices=masked_indices)
-        return get_total_spectrum_imzml(msi_data, num_bins=num_bins, masked_indices=masked_indices)
-    else # MzMLSource
-        return get_total_spectrum_mzml(msi_data, num_bins=num_bins, masked_indices=masked_indices)
-    end
 end
 
 """
@@ -1726,8 +1435,9 @@ function _iterate_compressed_fast(f::Function, data::MSIData, source::ImzMLSourc
         # For compressed data, we currently allocate new arrays per spectrum.
         # To truly minimize allocations, we'd need read_compressed_array! (in-place version).
         # For now, let's ensure we are at least using the optimized type-stable version.
-        mz_array = read_compressed_array(data, source.ibd_handle, meta.mz_asset, source.mz_format)
+        handle = get_handle(source)
-        intensity_array = read_compressed_array(data, source.ibd_handle, meta.int_asset, source.intensity_format)
+        mz_array = read_compressed_array(data, handle, meta.mz_asset, source.mz_format)
+        intensity_array = read_compressed_array(data, handle, meta.int_asset, source.intensity_format)
         
         mz_array .= ltoh.(mz_array)
         intensity_array .= ltoh.(intensity_array)
@@ -2044,19 +1754,19 @@ it will be created on-demand and cached.
 - A `BloomFilter{Int}` for the specified spectrum.
 """
 function get_bloom_filter(data::MSIData, index::Int)::BloomFilter{Int}
-    lock(data.cache_lock) do
+    # Double-checked locking for resilience if not precomputed
-        if data.bloom_filters === nothing
+    if data.bloom_filters === nothing || data.bloom_filters[index] === nothing
-            # Initialize Bloom filters on-demand
+        lock(data.cache_lock) do
-            data.bloom_filters = Vector{Union{BloomFilter{Int}, Nothing}}(undef, length(data.spectra_metadata))
+            if data.bloom_filters === nothing
-            fill!(data.bloom_filters, nothing) # Initialize with nothing
+                data.bloom_filters = Vector{Union{BloomFilter{Int}, Nothing}}(undef, length(data.spectra_metadata))
-        end
+                fill!(data.bloom_filters, nothing)
+            end
             
-        if data.bloom_filters[index] === nothing
+            if data.bloom_filters[index] === nothing
-            # Create Bloom filter lazily
+                mz, _ = GetSpectrum(data, index)
-            mz, intensity = GetSpectrum(data, index)
+                data.bloom_filters[index] = create_bloom_filter_for_spectrum(mz)
-            data.bloom_filters[index] = create_bloom_filter_for_spectrum(mz)
+            end
         end
-        
-        return data.bloom_filters[index]
     end
+    return data.bloom_filters[index]
 end

src/MSI_src.jl

@@ -21,7 +21,17 @@ export OpenMSIData,
        _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()
@@ -82,7 +92,6 @@ include("mzML.jl")
 include("imzML.jl")
 include("MzmlConverter.jl")
 include("Preprocessing.jl")
-include("FusedPipeline.jl")
 include("ImageProcessing.jl")
 include("Precalculations.jl")
 include("PreprocessingPipeline.jl")
@@ -111,13 +120,17 @@ This is the main entry point for the new data access API.
 - An `MSIData` object.
 """
 function OpenMSIData(filepath::String; cache_size=300, spectrum_type_map::Union{Dict{Int, Symbol}, Nothing}=nothing)
+    # Apply standard path normalization for cross-platform compatibility
+    # Ensure Windows backslashes are converted to OS-native separators
+    norm_filepath = normpath(replace(filepath, "\\" => "/"))
+    
     local msi_data
-    if endswith(lowercase(filepath), ".mzml")
+    if endswith(lowercase(norm_filepath), ".mzml")
-        msi_data = load_mzml_lazy(filepath, cache_size=cache_size)
+        msi_data = load_mzml_lazy(norm_filepath, cache_size=cache_size)
-    elseif endswith(lowercase(filepath), ".imzml")
+    elseif endswith(lowercase(norm_filepath), ".imzml")
-        msi_data = load_imzml_lazy(filepath, cache_size=cache_size)
+        msi_data = load_imzml_lazy(norm_filepath, cache_size=cache_size)
     else
-        error("Unsupported file type: $filepath. Please provide a .mzML or .imzML file.")
+        error("Unsupported file type: $norm_filepath. Please provide a .mzML or .imzML file.")
     end
 
     # Apply spectrum type map if provided

src/ResourcePool.jl

@@ -21,23 +21,7 @@ mutable struct ResourcePool{T}
     constructor::Function
 end
 
-"""
-    aligned_vector(::Type{T}, n::Int; alignment::Int=64) where T
 
-Creates a `Vector{T}` that is aligned to `alignment` bytes.
-Note: In modern Julia, standard vectors are often 16 or 64 byte aligned, but for
-HPC we ensure this by allocating slightly more and using a view, or using
-specific pointers. For simplicity and performance, we use a small hack:
-allocating a larger array and taking a 64-byte aligned view.
-"""
-function aligned_vector(::Type{T}, n::Int; alignment::Int=64) where T
-    # Allocate enough space to find an aligned starting point
-    raw = Vector{UInt8}(undef, n * sizeof(T) + alignment)
-    ptr = Int(pointer(raw))
-    off = (alignment - (ptr % alignment)) % alignment
-    # Return a reinterpret view of the aligned segment
-    return reinterpret(T, view(raw, (off + 1):(off + n * sizeof(T))))
-end
 
 """
     ResourcePool{T}(constructor::Function; max_size::Int=2 * nthreads())

src/imzML.jl

@@ -194,86 +194,6 @@ function axes_config_img(stream::IO)
     return param_groups
 end
 
-"""
-    get_spectrum_tag_offset(stream)
-
-Calculates the character offset within a `<spectrum>` tag, ignoring attribute values.
-
-# Arguments
-
-- `stream::IO`: The stream to parse.
-
-# Returns
-
-- `offset::Int`: The character offset.
-"""
-function get_spectrum_tag_offset(stream::IO)
-    offset = position(stream)
-    tag = find_tag(stream, r"^\s*<spectrum (.+)")
-    first = 1
-
-    while true
-        value = match(r"[^=]+\"([^\"]+)\"", tag.captures[1][first:end])
-        if value === nothing
-            break
-        end
-        first += value.offsets[1] + length(value.captures[1])
-        offset += length(value.captures[1])
-    end
-    return offset
-end
-
-"""
-    get_spectrum_attributes(stream, hIbd)
-
-Reads metadata to determine the byte offsets and data types for reading spectra.
-
-# Arguments
-
-- `stream::IO`: The stream to parse.
-- `hIbd::IO`: The IBD file handle.
-
-# Returns
-
-- `mz_first::Int`: The index of the first array
-- `array_lengths::Vector{Int}`: The lengths of the arrays.
-- `offset_mz::Int`: The byte offset of the mz array.
-- `offset_intensity::Int`: The byte offset of the intensity array.
-- `data_type_mz::DataType`: The data type of the mz array.
-- `data_type_intensity::DataType`: The data type of the intensity array.
-"""
-function get_spectrum_attributes(stream::IO, hIbd::IO)
-    # Look for position x
-    readuntil(stream, "IMS:1000050")
-    x_skip = 0
-    
-    # Look for position y  
-    readuntil(stream, "IMS:1000051")
-    y_skip = 0
-    
-    # Determine order of mz/intensity arrays
-    pos_before = position(stream)
-    
-    # Look for first external offset (could be mz or intensity)
-    readuntil(stream, "external offset")
-    current_line = readline(stream)
-    
-    # Check which array comes first by looking at the param group reference
-    mz_first = occursin("mzArray", current_line) ? 3 : 4
-    
-    # Find array length - look for the first external array length after coordinates
-    seek(stream, pos_before)
-    readuntil(stream, "IMS:1000103")
-    array_len_skip = 0
-    
-    # Find spectrum end
-    readuntil(stream, "</spectrum>")
-    spectrum_end_skip = 0
-    
-    return [x_skip, y_skip, mz_first, array_len_skip, spectrum_end_skip]
-end
-
-
 """
     read_spectrum_block(stream::IO)
 
@@ -454,7 +374,7 @@ end
 # Arguments
 
 - `stream::IO`: The stream to parse.
-- `hIbd::Union{IO, ThreadSafeFileHandle}`: The IBD file handle.
+- `hIbd::IO`: The IBD file handle.
 - `param_groups::Dict{String, SpecDim}`: The parameter groups.
 - `width::Int32`: The width of the image.
 - `height::Int32`: The height of the image.
@@ -469,7 +389,7 @@ end
 
 - `spectra_metadata::Vector{SpectrumMetadata}`: The spectrum metadata.
 """
-function parse_imzml_spectrum_block(stream::IO, hIbd::Union{IO, ThreadSafeFileHandle}, param_groups::Dict{String, SpecDim},
+function parse_imzml_spectrum_block(stream::IO, hIbd::IO, param_groups::Dict{String, SpecDim},
                                     width::Int32, height::Int32, num_spectra::Int32,
                                     default_mz_format::DataType, default_intensity_format::DataType, 
                                     mz_is_compressed::Bool, int_is_compressed::Bool, global_mode::SpectrumMode)
@@ -685,175 +605,6 @@ function load_imzml_lazy(file_path::String; cache_size::Int=100, use_mmap::Bool=
     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, 0.0, 0.0)
-            int_asset = SpectrumAsset(default_intensity_format, int_is_compressed, Int64(0), 0, :intensity, 0.0, 0.0)
-        else
-            mz_info = mz_data[1]
-            int_info = int_data[1]
-
-            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, 0.0, 0.0)
-            int_asset = SpectrumAsset(default_intensity_format, int_is_compressed, int_info.offset,
-                                     int_is_compressed ? int_info.encoded_length : int_info.array_length, :intensity, 0.0, 0.0)
-        end
-        
-        spectra_metadata[k] = SpectrumMetadata(x, y, "", :sample, spectrum_mode, mz_asset, int_asset)
-    end
-
-    return spectra_metadata
-end
-
 
 # =============================================================================
 #

test/runtests.jl

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