Julian/JuliaMSI
1
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

Added bloom filters custom API, thread safety for multithreading several important functions, streaming processing to reduce GC and Memory allocations calls, buffer pool changed with lock safety, and debug memory monitoring to track leaks, alongside better loading times in general

Browse Source
This commit is contained in:
Pixelguy14 2025-11-09 13:41:29 -06:00
parent a8c2dfcc91
commit edfbb77f58
11 changed files with 1256 additions and 491 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

12
Manifest.toml
View File

@ -2,7 +2,7 @@
julia_version = "1.11.7"
manifest_format = "2.0"
project_hash = "f5ab935406d0b7b7e657445ecfa45944ea71f247"
project_hash = "38f13d997585a9af185f1863ff626d584a7e791e"
[[deps.ATK_jll]]
deps = ["Artifacts", "Glib_jll", "JLLWrappers", "Libdl"]
@ -651,6 +651,12 @@ git-tree-sha1 = "05882d6995ae5c12bb5f36dd2ed3f61c98cbb172"
uuid = "53c48c17-4a7d-5ca2-90c5-79b7896eea93"
version = "0.8.5"
[[deps.FlameGraphs]]
deps = ["AbstractTrees", "Colors", "FileIO", "FixedPointNumbers", "IndirectArrays", "LeftChildRightSiblingTrees", "Profile"]
git-tree-sha1 = "0166baf81babb91cf78bfcc771d8e87c43d568df"
uuid = "08572546-2f56-4bcf-ba4e-bab62c3a3f89"
version = "1.1.0"
[[deps.Fontconfig_jll]]
deps = ["Artifacts", "Bzip2_jll", "Expat_jll", "FreeType2_jll", "JLLWrappers", "Libdl", "Libuuid_jll", "Zlib_jll"]
git-tree-sha1 = "f85dac9a96a01087df6e3a749840015a0ca3817d"
@ -1779,6 +1785,10 @@ deps = ["Unicode"]
uuid = "de0858da-6303-5e67-8744-51eddeeeb8d7"
version = "1.11.0"
[[deps.Profile]]
uuid = "9abbd945-dff8-562f-b5e8-e1ebf5ef1b79"
version = "1.11.0"
[[deps.ProgressMeter]]
deps = ["Distributed", "Printf"]
git-tree-sha1 = "fbb92c6c56b34e1a2c4c36058f68f332bec840e7"

42
Project.toml
View File

@ -12,6 +12,7 @@ Colors = "5ae59095-9a9b-59fe-a467-6f913c188581"
DataFrames = "a93c6f00-e57d-5684-b7b6-d8193f3e46c0"
Dates = "ade2ca70-3891-5945-98fb-dc099432e06a"
FileIO = "5789e2e9-d7fb-5bc7-8068-2c6fae9b9549"
FlameGraphs = "08572546-2f56-4bcf-ba4e-bab62c3a3f89"
GLMakie = "e9467ef8-e4e7-5192-8a1a-b1aee30e663a"
Genie = "c43c736e-a2d1-11e8-161f-af95117fbd1e"
GenieFramework = "a59fdf5c-6bf0-4f5d-949c-a137c9e2f353"
@ -19,6 +20,7 @@ HistogramThresholding = "2c695a8d-9458-5d45-9878-1b8a99cf7853"
ImageBinarization = "cbc4b850-ae4b-5111-9e64-df94c024a13d"
ImageComponentAnalysis = "d9b9e9a0-1569-11e9-2cb5-bbca914b0e89"
ImageContrastAdjustment = "f332f351-ec65-5f6a-b3d1-319c6670881a"
ImageCore = "a09fc81d-aa75-5fe9-8630-4744c3626534"
ImageFiltering = "6a3955dd-da59-5b1f-98d4-e7296123deb5"
ImageMorphology = "787d08f9-d448-5407-9aad-5290dd7ab264"
ImageSegmentation = "80713f31-8817-5129-9cf8-209ff8fb23e1"
@ -39,3 +41,43 @@ Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"
StatsBase = "2913bbd2-ae8a-5f71-8c99-4fb6c76f3a91"
StipplePlotly = "ec984513-233d-481d-95b0-a3b58b97af2b"
UUIDs = "cf7118a7-6976-5b1a-9a39-7adc72f591a4"
[compat]
Accessors = "0.1"
Base64 = "1.11"
CSV = "0.10"
CairoMakie = "0.13"
ColorSchemes = "3.31"
Colors = "0.12"
DataFrames = "1.7"
Dates = "1.11"
FileIO = "1.17"
GLMakie = "0.11"
Genie = "5.31"
GenieFramework = "2.8"
HistogramThresholding = "0.3"
ImageBinarization = "0.3"
ImageComponentAnalysis = "0.2"
ImageContrastAdjustment = "0.3"
ImageCore = "0.10.5"
ImageFiltering = "0.7"
ImageMorphology = "0.4"
ImageSegmentation = "1.9"
Images = "0.26"
JSON = "0.21"
Libz = "1.0"
LinearAlgebra = "1.11"
Loess = "0.6"
Mmap = "1.11"
NativeFileDialog = "0.2"
NaturalSort = "1.0"
PlotlyBase = "0.8"
Printf = "1.11"
ProgressMeter = "1.11"
SavitzkyGolay = "0.9"
Serialization = "1.11"
Statistics = "1.11"
StatsBase = "0.34"
StipplePlotly = "0.13"
UUIDs = "1.11"
julia = "1.11"

85
app.jl
View File

@ -19,9 +19,10 @@ using Base.Filesystem: mv # To rename files in the system
using Printf # Required for @sprintf macro in colorbar generation
using JSON
using Dates
using Base.Threads
# Bring MSIData into App module's scope
using .MSI_src: MSIData, OpenMSIData, process_spectrum, IterateSpectra, ImzMLSource, _iterate_spectra_fast, MzMLSource, find_mass, ViridisPalette, get_mz_slice, get_multiple_mz_slices, quantize_intensity, save_bitmap, median_filter, save_bitmap, downsample_spectrum, TrIQ, precompute_analytics, ImportMzmlFile, generate_colorbar_image, load_and_prepare_mask
using .MSI_src: MSIData, OpenMSIData, process_spectrum, IterateSpectra, ImzMLSource, _iterate_spectra_fast, MzMLSource, find_mass, ViridisPalette, get_mz_slice, get_multiple_mz_slices, quantize_intensity, save_bitmap, median_filter, save_bitmap, downsample_spectrum, TrIQ, precompute_analytics, ImportMzmlFile, generate_colorbar_image, load_and_prepare_mask, set_global_mz_range!
if !@isdefined(increment_image)
include("./julia_imzML_visual.jl")
@ -837,7 +838,7 @@ end
end
end
@onbutton createMeanPlot begin
@onbutton createMeanPlot @time begin
if isempty(selected_folder_main)
msg = "No dataset selected. Please process a file and select a folder first."
warning_msg = true
@ -862,12 +863,18 @@ end
end
if msi_data === nothing || full_route != target_path
if msi_data !== nothing
close(msi_data)
end
msg = "Reloading $(basename(target_path)) for analysis..."
full_route = target_path
msi_data = OpenMSIData(target_path)
if haskey(get(entry, "metadata", Dict()), "global_min_mz") && entry["metadata"]["global_min_mz"] !== nothing
msi_data.global_min_mz = entry["metadata"]["global_min_mz"]
msi_data.global_max_mz = entry["metadata"]["global_max_mz"]
raw_min = entry["metadata"]["global_min_mz"]
raw_max = entry["metadata"]["global_max_mz"]
min_val = isa(raw_min, Dict) ? get(raw_min, "value", raw_min) : raw_min
max_val = isa(raw_max, Dict) ? get(raw_max, "value", raw_max) : raw_max
set_global_mz_range!(msi_data, convert(Float64, min_val), convert(Float64, max_val))
else
precompute_analytics(msi_data)
end
@ -904,7 +911,7 @@ end
end
end
@onbutton createSumPlot begin
@onbutton createSumPlot @time begin
if isempty(selected_folder_main)
msg = "No dataset selected. Please process a file and select a folder first."
warning_msg = true
@ -914,32 +921,37 @@ end
progressSpectraPlot = true
btnPlotDisable = true
btnStartDisable = true
msg = "Loading total spectrum plot for $(selected_folder_main)..."
msg = "Loading total spectrum plot for $(selected_folder_main)..."
try
sTime = time()
registry = load_registry(registry_path)
entry = registry[selected_folder_main]
target_path = entry["source_path"]
try
sTime = time()
registry = load_registry(registry_path)
entry = registry[selected_folder_main]
target_path = entry["source_path"]
if target_path == "unknown (manually added)"
msg = "Dataset selected contained no route."
warning_msg = true
return
end
if msi_data === nothing || full_route != target_path
msg = "Reloading $(basename(target_path)) for analysis..."
full_route = target_path
msi_data = OpenMSIData(target_path)
if haskey(get(entry, "metadata", Dict()), "global_min_mz") && entry["metadata"]["global_min_mz"] !== nothing
msi_data.global_min_mz = entry["metadata"]["global_min_mz"]
msi_data.global_max_mz = entry["metadata"]["global_max_mz"]
else
precompute_analytics(msi_data)
end
end
if target_path == "unknown (manually added)"
msg = "Dataset selected contained no route."
warning_msg = true
return
end
if msi_data === nothing || full_route != target_path
if msi_data !== nothing
close(msi_data)
end
msg = "Reloading $(basename(target_path)) for analysis..."
full_route = 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!(msi_data, convert(Float64, min_val), convert(Float64, max_val))
else
precompute_analytics(msi_data)
end
end
local mask_path_for_plot::Union{String, Nothing} = nothing
if maskEnabled && get(entry, "has_mask", false)
mask_path_for_plot = get(entry, "mask_path", "")
@ -971,7 +983,7 @@ end
end
end
@onbutton createXYPlot begin
@onbutton createXYPlot @time begin
if isempty(selected_folder_main)
msg = "No dataset selected. Please process a file and select a folder first."
warning_msg = true
@ -1005,12 +1017,18 @@ end
end
if msi_data === nothing || full_route != target_path
if msi_data !== nothing
close(msi_data)
end
msg = "Reloading $(basename(target_path)) for analysis..."
full_route = target_path
msi_data = OpenMSIData(target_path)
if haskey(get(entry, "metadata", Dict()), "global_min_mz") && entry["metadata"]["global_min_mz"] !== nothing
msi_data.global_min_mz = entry["metadata"]["global_min_mz"]
msi_data.global_max_mz = entry["metadata"]["global_max_mz"]
raw_min = entry["metadata"]["global_min_mz"]
raw_max = entry["metadata"]["global_max_mz"]
min_val = isa(raw_min, Dict) ? get(raw_min, "value", raw_min) : raw_min
max_val = isa(raw_max, Dict) ? get(raw_max, "value", raw_max) : raw_max
set_global_mz_range!(msi_data, convert(Float64, min_val), convert(Float64, max_val))
else
precompute_analytics(msi_data)
end
@ -1351,6 +1369,9 @@ end
# This handler will now correctly load the first image from the newly selected folder.
@onchange selected_folder_main begin
if msi_data !== nothing
close(msi_data)
end
msi_data = nothing
log_memory_usage("Folder Changed (msi_data cleared)", msi_data)
GC.gc()
@ -1900,7 +1921,7 @@ end
@mounted watchplots()
@onchange isready begin
@onchange isready @time begin
if isready && !registry_init_done
warmup_init()
try

4
julia_imzML_visual.jl
View File

@ -847,8 +847,8 @@ function extract_metadata(msi_data::MSIData, source_path::String)
Dict("parameter" => "File Name", "value" => basename(source_path)),
Dict("parameter" => "Number of Spectra", "value" => length(msi_data.spectra_metadata)),
Dict("parameter" => "Image Dimensions", "value" => "$(msi_data.image_dims[1]) x $(msi_data.image_dims[2])"),
Dict("parameter" => "Global Min m/z", "value" => @sprintf("%.4f", msi_data.global_min_mz)),
Dict("parameter" => "Global Max m/z", "value" => @sprintf("%.4f", msi_data.global_max_mz)),
Dict("parameter" => "Global Min m/z", "value" => @sprintf("%.4f", Threads.atomic_add!(msi_data.global_min_mz, 0.0))),
Dict("parameter" => "Global Max m/z", "value" => @sprintf("%.4f", Threads.atomic_add!(msi_data.global_max_mz, 0.0))),
Dict("parameter" => "Mean TIC", "value" => @sprintf("%.2e", mean(df.TIC))),
Dict("parameter" => "Mean BPI", "value" => @sprintf("%.2e", mean(df.BPI))),
Dict("parameter" => "Mean # Points", "value" => @sprintf("%.1f", mean(df.NumPoints))),

205
src/BloomFilters.jl Normal file
View File

@ -0,0 +1,205 @@
# src/BloomFilter.jl
"""
BloomFilter{T}
A simple, efficient Bloom filter implementation for probabilistic set membership testing.
# 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)
"""
mutable struct BloomFilter{T}
bits::BitVector
size::Int
hash_count::Int
seed::UInt64
count::Int
end
"""
StreamingBloomFilter
A memory-efficient Bloom filter that doesn't store all values at once.
"""
mutable struct StreamingBloomFilter
bits::BitVector
size::Int
hash_count::Int
seed::UInt64
count::Int
# Streaming state
current_chunk::Vector{Int}
chunk_size::Int
end
"""
BloomFilter(expected_elements::Int, false_positive_rate::Float64=0.01; seed::UInt64=0x12345678)
Creates a Bloom filter optimized for the expected number of elements and desired false positive rate.
"""
function BloomFilter{T}(expected_elements::Int, false_positive_rate::Float64=0.01; seed::Union{UInt32,UInt64}=0x12345678) where T
# Convert seed to UInt64 for consistency
seed_uint64 = UInt64(seed)
# Calculate optimal parameters using standard formulas
size = optimal_bit_size(expected_elements, false_positive_rate)
hash_count = optimal_hash_count(expected_elements, size)
bits = falses(size)
return BloomFilter{T}(bits, size, hash_count, seed_uint64, 0)
end
"""
optimal_bit_size(n::Int, p::Float64) -> Int
Calculates the optimal number of bits for a Bloom filter given:
- n: expected number of elements
- p: desired false positive rate
"""
function optimal_bit_size(n::Int, p::Float64)::Int
if p <= 0.0 || p >= 1.0
throw(ArgumentError("False positive rate must be between 0 and 1"))
end
# m = - (n * ln(p)) / (ln(2)^2)
m = ceil(Int, - (n * log(p)) / (log(2)^2))
return max(m, 1)
end
"""
optimal_hash_count(n::Int, m::Int) -> Int
Calculates the optimal number of hash functions for a Bloom filter.
"""
function optimal_hash_count(n::Int, m::Int)::Int
if n <= 0 || m <= 0
return 1
end
# k = (m / n) * ln(2)
k = max(1, round(Int, (m / n) * log(2)))
return min(k, 8) # Practical limit to avoid too many hashes
end
"""
hash_functions(item::T, count::Int, size::Int, seed::UInt64) -> Vector{Int}
Generates multiple hash values for an item using double hashing technique.
"""
function hash_functions(item::T, count::Int, size::Int, seed::UInt64) where T
# Use Julia's built-in hash with different seeds
hashes = Vector{Int}(undef, count)
# First hash with the main seed
h1 = hash(item, seed)
h2 = hash(item, seed + 1)
for i in 1:count
# Double hashing: h_i = h1 + i * h2
combined_hash = UInt64(h1) + UInt64(i) * UInt64(h2)
hashes[i] = (combined_hash % UInt64(size)) + 1 # 1-based indexing
end
return hashes
end
"""
Base.push!(bf::BloomFilter{T}, item::T)
Adds an element to the Bloom filter.
"""
function Base.push!(bf::BloomFilter{T}, item::T) where T
hashes = hash_functions(item, bf.hash_count, bf.size, bf.seed)
for h in hashes
bf.bits[h] = true
end
bf.count += 1
return bf
end
"""
Base.in(item::T, bf::BloomFilter{T}) -> Bool
Checks whether an element is possibly in the Bloom filter.
Returns `true` if the element might be in the set, `false` if it's definitely not.
"""
function Base.in(item::T, bf::BloomFilter{T})::Bool where T
hashes = hash_functions(item, bf.hash_count, bf.size, bf.seed)
for h in hashes
if !bf.bits[h]
return false # Definitely not in the set
end
end
return true # Possibly in the set
end
"""
contains(bf::BloomFilter{T}, item::T) -> Bool
Alias for `in()` for compatibility with your existing code.
"""
contains(bf::BloomFilter{T}, item::T) where T = item in bf
"""
add!(bf::BloomFilter{T}, item::T)
Alias for `push!()` for compatibility with your existing code.
"""
add!(bf::BloomFilter{T}, item::T) where T = push!(bf, item)
"""
false_positive_rate(bf::BloomFilter) -> Float64
Estimates the current false positive rate of the Bloom filter.
"""
function false_positive_rate(bf::BloomFilter)::Float64
if bf.count == 0
return 0.0
end
# Theoretical false positive rate: (1 - e^(-k * n / m)) ^ k
k = bf.hash_count
n = bf.count
m = bf.size
return (1 - exp(-k * n / m)) ^ k
end
"""
fill_ratio(bf::BloomFilter) -> Float64
Returns the fraction of bits that are set to 1.
"""
function fill_ratio(bf::BloomFilter)::Float64
return count(bf.bits) / length(bf.bits)
end
"""
is_empty(bf::BloomFilter) -> Bool
Checks if the Bloom filter is empty (no elements added).
"""
function is_empty(bf::BloomFilter)::Bool
return bf.count == 0
end
"""
reset!(bf::BloomFilter)
Clears the Bloom filter, removing all elements.
"""
function reset!(bf::BloomFilter)
fill!(bf.bits, false)
bf.count = 0
return bf
end
# Specialized constructor for empty Bloom filters
function BloomFilter{T}(;size::Int=100, hash_count::Int=3, seed::Union{UInt32,UInt64}=0x12345678) where T
seed_uint64 = UInt64(seed)
bits = falses(size)
return BloomFilter{T}(bits, size, hash_count, seed_uint64, 0)
end

185
src/Common.jl
View File

@ -1,4 +1,114 @@
# src/Common.jl
# src/Common.jl - Updated with BloomFilter
using Base.Threads
# --- Buffer Pooling ---
"""
BufferPool
A thread-safe pool of `Vector{UInt8}` buffers, categorized by their size.
This helps reduce memory allocations and garbage collection overhead by reusing buffers.
# Fields
- `lock`: A `ReentrantLock` to ensure thread-safe access to the pool.
- `buffers`: A dictionary mapping buffer size (in bytes) to a list of available buffers of that size.
"""
mutable struct BufferPool
lock::ReentrantLock
buffers::Dict{Int, Vector{Vector{UInt8}}}
function BufferPool()
new(ReentrantLock(), Dict{Int, Vector{Vector{UInt8}}}())
end
end
"""
get_buffer!(pool::BufferPool, size::Int) -> Vector{UInt8}
Retrieves a `Vector{UInt8}` buffer of at least `size` bytes from the pool.
If no suitable buffer is available, a new one is allocated.
# Arguments
- `pool`: The `BufferPool` to retrieve the buffer from.
- `size`: The minimum desired size of the buffer in bytes.
# Returns
- A `Vector{UInt8}` buffer.
"""
function get_buffer!(pool::BufferPool, size::Int)::Vector{UInt8}
lock(pool.lock) do
# Find an existing buffer that is large enough
for (buffer_size, buffer_list) in pool.buffers
if buffer_size >= size && !isempty(buffer_list)
buffer = pop!(buffer_list)
# Resize if necessary (and if it's significantly larger than needed)
if length(buffer) > size * 2 # Heuristic: if buffer is more than twice the requested size
return Vector{UInt8}(undef, size) # Allocate new smaller buffer
else
return buffer
end
end
end
# No suitable buffer found, allocate a new one
return Vector{UInt8}(undef, size)
end
end
"""
release_buffer!(pool::BufferPool, buffer::Vector{UInt8})
Returns a `Vector{UInt8}` buffer to the pool for future reuse.
# Arguments
- `pool`: The `BufferPool` to return the buffer to.
- `buffer`: The `Vector{UInt8}` buffer to release.
"""
function release_buffer!(pool::BufferPool, buffer::Vector{UInt8})
lock(pool.lock) do
buffer_size = length(buffer)
if !haskey(pool.buffers, buffer_size)
pool.buffers[buffer_size] = Vector{Vector{UInt8}}()
end
push!(pool.buffers[buffer_size], buffer)
end
return nothing
end
"""
SimpleBufferPool
A dramatically simplified buffer pool that avoids complex locking.
"""
mutable struct SimpleBufferPool
buffers::Dict{Int, Vector{Vector{UInt8}}}
max_pool_size::Int
end
SimpleBufferPool() = SimpleBufferPool(Dict{Int, Vector{Vector{UInt8}}}(), 50)
function get_buffer!(pool::SimpleBufferPool, size::Int)::Vector{UInt8}
# Check for existing buffers of exact size first
if haskey(pool.buffers, size) && !isempty(pool.buffers[size])
return pop!(pool.buffers[size])
end
# No suitable buffer found, allocate new one
return Vector{UInt8}(undef, size)
end
function release_buffer!(pool::SimpleBufferPool, buffer::Vector{UInt8})
size = length(buffer)
if !haskey(pool.buffers, size)
pool.buffers[size] = Vector{Vector{UInt8}}()
end
# Limit pool size to prevent memory bloat
if length(pool.buffers[size]) < pool.max_pool_size
push!(pool.buffers[size], buffer)
end
# If pool is full, let buffer get GC'd
end
# --- Unified Error Types ---
abstract type MSIError <: Exception end
@ -19,19 +129,13 @@ end
# --- Shared Constants ---
const DEFAULT_CACHE_SIZE = 100
const DEFAULT_NUM_BINS = 2000
const DEFAULT_BLOOM_FILTER_SIZE = 1000 # Default size for empty spectra
const DEFAULT_FALSE_POSITIVE_RATE = 0.01
"""
validate_spectrum_data(mz, intensity, id)
Checks a spectrum for common data integrity issues.
- Checks for `NaN` or `Inf` values in intensity and m/z arrays.
- Verifies that the m/z array is sorted in ascending order.
- Ensures that m/z and intensity arrays have the same length.
# Throws
- `InvalidSpectrumError` if any of the checks fail.
"""
function validate_spectrum_data(mz::AbstractVector, intensity::AbstractVector, id)
if length(mz) != length(intensity)
@ -48,3 +152,66 @@ function validate_spectrum_data(mz::AbstractVector, intensity::AbstractVector, i
return true
end
"""
create_bloom_filter_for_spectrum(mz::AbstractVector{<:Real}) -> BloomFilter{Int}
Memory-optimized version that processes data in chunks to reduce temporary allocations
and avoid holding large intermediate arrays.
"""
function create_bloom_filter_for_spectrum(mz::AbstractVector{<:Real})::BloomFilter{Int}
if isempty(mz)
return empty_bloom_filter()
end
# Size Bloom filter appropriately
n_points = length(mz)
bf = BloomFilter{Int}(n_points, DEFAULT_FALSE_POSITIVE_RATE)
discretization_factor = 100.0
chunk_size = 10000 # Process in chunks to reduce memory pressure
# Process spectrum data in chunks to avoid large temporary arrays
for chunk_start in 1:chunk_size:n_points
chunk_end = min(chunk_start + chunk_size - 1, n_points)
@inbounds for i in chunk_start:chunk_end
discretized_val = round(Int, mz[i] * discretization_factor)
push!(bf, discretized_val)
end
end
return bf
end
"""
empty_bloom_filter() -> BloomFilter{Float64}
Creates an empty Bloom filter for spectra with no data.
"""
function empty_bloom_filter()::BloomFilter{Int}
return BloomFilter{Int}(size=DEFAULT_BLOOM_FILTER_SIZE, hash_count=3)
end
"""
AtomicFlag
A thread-safe boolean flag using atomic operations.
"""
struct AtomicFlag
value::Base.Threads.Atomic{Int}
end
AtomicFlag() = AtomicFlag(Base.Threads.Atomic{Int}(0))
function set!(flag::AtomicFlag)
Base.Threads.atomic_add!(flag.value, 1)
end
function is_set(flag::AtomicFlag)::Bool
return Base.Threads.atomic_add!(flag.value, 0) > 0
end
function reset!(flag::AtomicFlag)
Base.Threads.atomic_xchg!(flag.value, 0)
end

658
src/MSIData.jl
View File

@ -8,6 +8,64 @@ efficiently.
using Base64, Libz, Serialization, Printf, DataFrames, Base.Threads # For reading binary data
const FILE_HANDLE_LOCK = ReentrantLock()
"""
ThreadSafeFileHandle
Wraps a file handle with thread-safe operations.
"""
mutable struct ThreadSafeFileHandle
path::String
handle::IO
lock::ReentrantLock
end
function ThreadSafeFileHandle(path::String, mode::String="r")
handle = lock(FILE_HANDLE_LOCK) do
open(path, mode)
end
return ThreadSafeFileHandle(path, handle, ReentrantLock())
end
function Base.seek(tsfh::ThreadSafeFileHandle, pos::Integer)
lock(tsfh.lock) do
seek(tsfh.handle, pos)
end
end
function Base.read(tsfh::ThreadSafeFileHandle, n::Integer)
lock(tsfh.lock) do
return read(tsfh.handle, n)
end
end
function Base.read!(tsfh::ThreadSafeFileHandle, a::AbstractArray)
lock(tsfh.lock) do
return read!(tsfh.handle, a)
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
# Abstract type for different data sources (e.g., mzML, imzML)
# This allows dispatching to the correct binary reading logic.
abstract type MSDataSource end
@ -20,7 +78,7 @@ A data source for `.imzML` files, holding a handle to the binary `.ibd` file
and the expected format for m/z and intensity arrays.
"""
struct ImzMLSource <: MSDataSource
ibd_handle::IO
ibd_handle::Union{IO, ThreadSafeFileHandle}
mz_format::Type
intensity_format::Type
end
@ -32,7 +90,7 @@ A data source for `.mzML` files, holding a handle to the `.mzML` file itself
(which contains the binary data encoded in Base64) and the expected data formats.
"""
struct MzMLSource <: MSDataSource
file_handle::IO
file_handle::Union{IO, ThreadSafeFileHandle}
mz_format::Type
intensity_format::Type
end
@ -122,24 +180,31 @@ efficient repeated access to spectra.
mutable struct MSIData
source::MSDataSource
spectra_metadata::Vector{SpectrumMetadata}
image_dims::Tuple{Int, Int} # (width, height) for imaging data
coordinate_map::Union{Matrix{Int}, Nothing} # Maps (x,y) to linear index for imzML
image_dims::Tuple{Int, Int}
coordinate_map::Union{Matrix{Int}, Nothing}
# LRU Cache for GetSpectrum
cache::Dict{Int, Tuple{Vector{Float64}, Vector{Float64}}}
cache_order::Vector{Int} # Stores indices, with most recently used at the end
cache_size::Int # Max number of spectra in cache
cache_lock::ReentrantLock # To make cache access thread-safe
cache_order::Vector{Int}
cache_size::Int
cache_lock::ReentrantLock
# Pre-computed analytics/metadata
global_min_mz::Union{Float64, Nothing}
global_max_mz::Union{Float64, Nothing}
# Buffer Pool for binary data operations
buffer_pool::SimpleBufferPool
# Pre-computed analytics/metadata - use Threads.Atomic for compatibility
global_min_mz::Threads.Atomic{Float64}
global_max_mz::Threads.Atomic{Float64}
spectrum_stats_df::Union{DataFrame, Nothing}
bloom_filters::Union{Vector{<:BloomFilter}, Nothing}
analytics_ready::AtomicFlag
function MSIData(source, metadata, dims, coordinate_map, cache_size)
obj = new(source, metadata, dims, coordinate_map,
Dict(), [], cache_size, ReentrantLock(),
nothing, nothing, nothing) # Initialize new fields to nothing
SimpleBufferPool(),
Threads.Atomic{Float64}(Inf), Threads.Atomic{Float64}(-Inf),
nothing, nothing, AtomicFlag())
# Ensure file handles are closed when the object is garbage collected
finalizer(obj) do o
@ -153,6 +218,38 @@ mutable struct MSIData
end
end
# Thread-safe state management
function set_global_mz_range!(data::MSIData, min_mz::Float64, max_mz::Float64)
Threads.atomic_xchg!(data.global_min_mz, min_mz)
Threads.atomic_xchg!(data.global_max_mz, max_mz)
end
function get_global_mz_range(data::MSIData)
# Atomic read using atomic_add! with zero
min_mz = Threads.atomic_add!(data.global_min_mz, 0.0)
max_mz = Threads.atomic_add!(data.global_max_mz, 0.0)
return (min_mz, max_mz)
end
function set_analytics_data!(data::MSIData, stats_df::DataFrame, bloom_filters::Vector)
lock(data.cache_lock) do
data.spectrum_stats_df = stats_df
data.bloom_filters = bloom_filters
end
end
function get_bloom_filters(data::MSIData)
lock(data.cache_lock) do
return data.bloom_filters
end
end
function get_spectrum_stats(data::MSIData)
lock(data.cache_lock) do
return data.spectrum_stats_df
end
end
"""
Base.close(data::MSIData)
@ -286,21 +383,39 @@ by this function and is assumed to be handled by the caller if necessary.
# Returns
- A `Vector` of the appropriate type containing the decoded data.
"""
function read_binary_vector(io::IO, asset::SpectrumAsset)
function read_binary_vector(data::MSIData, io::IO, asset::SpectrumAsset)
if asset.offset < 0 || asset.offset >= filesize(io)
throw(FileFormatError("Invalid asset offset: $(asset.offset) for file size $(filesize(io))"))
end
seek(io, asset.offset)
raw_b64 = read(io, asset.encoded_length)
decoded_bytes = base64decode(strip(String(raw_b64)))
bytes_io = IOBuffer(asset.is_compressed ? Libz.inflate(decoded_bytes) : decoded_bytes)
n_elements = Int(bytes_io.size / sizeof(asset.format))
out_array = Array{asset.format}(undef, n_elements)
read!(bytes_io, out_array)
# Use String directly to avoid intermediate allocations
b64_string = String(raw_b64)
# FIX: little-endian to host byte order for mzML data
out_array .= ltoh.(out_array)
local decoded_bytes::Vector{UInt8}
if asset.is_compressed
# Direct Base64 decode to temporary, then decompress
temp_decoded = Base64.base64decode(b64_string)
decoded_bytes = Libz.inflate(temp_decoded)
else
# Direct Base64 decode
decoded_bytes = Base64.base64decode(b64_string)
end
# Calculate number of elements
n_elements = length(decoded_bytes) ÷ sizeof(asset.format)
out_array = Vector{asset.format}(undef, n_elements)
# Use unsafe_wrap to avoid copy - this is the critical optimization
temp_array = unsafe_wrap(Vector{asset.format}, pointer(decoded_bytes), n_elements)
# Convert byte order in-place
@inbounds for i in 1:n_elements
out_array[i] = ltoh(temp_array[i])
end
return out_array
end
@ -344,11 +459,9 @@ function read_spectrum_from_disk(source::ImzMLSource, meta::SpectrumMetadata)
mz .= ltoh.(mz)
intensity .= ltoh.(intensity) # FIX: Added missing conversion for intensity
mz_f64, intensity_f64 = Float64.(mz), Float64.(intensity)
validate_spectrum_data(mz, intensity, meta.id)
validate_spectrum_data(mz_f64, intensity_f64, meta.id)
return mz_f64, intensity_f64
return mz, intensity
end
"""
@ -365,16 +478,16 @@ This is a high-level function that orchestrates the reading process by calling
# Returns
- A tuple `(mz, intensity)` containing the two requested data arrays.
"""
function read_spectrum_from_disk(source::MzMLSource, meta::SpectrumMetadata)
function read_spectrum_from_disk(data::MSIData, source::MzMLSource, meta::SpectrumMetadata)
# For mzML, data is Base64 encoded within the XML
mz = read_binary_vector(source.file_handle, meta.mz_asset)
intensity = read_binary_vector(source.file_handle, meta.int_asset)
mz = read_binary_vector(data, source.file_handle, meta.mz_asset)
intensity = read_binary_vector(data, source.file_handle, meta.int_asset)
mz_f64, intensity_f64 = Float64.(mz), Float64.(intensity)
validate_spectrum_data(mz_f64, intensity_f64, meta.id)
validate_spectrum_data(mz, intensity, meta.id)
return mz_f64, intensity_f64
return mz, intensity
end
# --- Public API --- #
@ -414,7 +527,13 @@ function GetSpectrum(data::MSIData, index::Int)
# Phase 2: Cache Miss - Read from disk (no lock)
meta = data.spectra_metadata[index]
spectrum = read_spectrum_from_disk(data.source, meta)
local spectrum
if data.source isa ImzMLSource
spectrum = read_spectrum_from_disk(data.source, meta)
else # MzMLSource
spectrum = read_spectrum_from_disk(data, data.source, meta)
end
# Phase 3: Update cache (with lock) and get final value
return lock(data.cache_lock) do
@ -527,122 +646,118 @@ as they can use this cached data. This function modifies the `MSIData` object in
and is idempotent.
"""
function precompute_analytics(msi_data::MSIData)
# Idempotency check
if msi_data.spectrum_stats_df !== nothing && hasproperty(msi_data.spectrum_stats_df, :MinMZ)
println("Analytics have already been pre-computed.")
# Double-checked locking pattern
if is_set(msi_data.analytics_ready)
return
end
println("Pre-computing analytics (single pass)...")
start_time = time_ns()
num_spectra = length(msi_data.spectra_metadata)
# DEBUG: Check the first spectrum's metadata
if num_spectra > 0
first_meta = msi_data.spectra_metadata[1]
println("DEBUG First spectrum metadata:")
println(" mz_asset: format=$(first_meta.mz_asset.format), compressed=$(first_meta.mz_asset.is_compressed)")
println(" mz_asset: offset=$(first_meta.mz_asset.offset), encoded_length=$(first_meta.mz_asset.encoded_length)")
println(" int_asset: format=$(first_meta.int_asset.format), compressed=$(first_meta.int_asset.is_compressed)")
println(" int_asset: offset=$(first_meta.int_asset.offset), encoded_length=$(first_meta.int_asset.encoded_length)")
println(" mode: $(first_meta.mode)")
end
# Initialize variables for global stats
g_min_mz = Inf
g_max_mz = -Inf
# Initialize vectors for per-spectrum stats
tics = Vector{Float64}(undef, num_spectra)
bpis = Vector{Float64}(undef, num_spectra)
bp_mzs = Vector{Float64}(undef, num_spectra)
num_points = Vector{Int}(undef, num_spectra)
min_mzs = Vector{Float64}(undef, num_spectra)
max_mzs = Vector{Float64}(undef, num_spectra)
modes = Vector{SpectrumMode}(undef, num_spectra)
is_compressed = Vector{Bool}(undef, num_spectra)
# DEBUG: Add counter to see how many spectra have data
spectra_with_data = 0
empty_spectra = 0
_iterate_spectra_fast(msi_data) do idx, mz, intensity
# Store metadata
modes[idx] = msi_data.spectra_metadata[idx].mode
is_compressed[idx] = msi_data.spectra_metadata[idx].mz_asset.is_compressed ||
msi_data.spectra_metadata[idx].int_asset.is_compressed
if isempty(mz)
empty_spectra += 1
tics[idx] = 0.0
bpis[idx] = 0.0
bp_mzs[idx] = 0.0
num_points[idx] = 0
min_mzs[idx] = Inf
max_mzs[idx] = -Inf
lock(msi_data.cache_lock) do
if is_set(msi_data.analytics_ready)
return
else
spectra_with_data += 1
end
# Update global m/z range
local_min, local_max = extrema(mz)
g_min_mz = min(g_min_mz, local_min)
g_max_mz = max(g_max_mz, local_max)
min_mzs[idx] = local_min
max_mzs[idx] = local_max
println("Pre-computing analytics (streaming mode)...")
start_time = time_ns()
# Calculate per-spectrum stats
tics[idx] = sum(intensity)
max_int, max_idx = findmax(intensity)
bpis[idx] = max_int
bp_mzs[idx] = mz[max_idx]
num_points[idx] = length(mz)
num_spectra = length(msi_data.spectra_metadata)
# Initialize thread-local variables for global stats
thread_local_min_mz = [Inf for _ in 1:Threads.nthreads()]
thread_local_max_mz = [-Inf for _ in 1:Threads.nthreads()]
# Initialize vectors for per-spectrum stats
tics = Vector{Float64}(undef, num_spectra)
bpis = Vector{Float64}(undef, num_spectra)
bp_mzs = Vector{Float64}(undef, num_spectra)
num_points = Vector{Int}(undef, num_spectra)
min_mzs = Vector{Float64}(undef, num_spectra)
max_mzs = Vector{Float64}(undef, num_spectra)
modes = Vector{SpectrumMode}(undef, num_spectra)
# Don't precompute Bloom filters by default - create them lazily
bloom_filters = nothing
# Process in chunks to reduce memory pressure
chunk_size = min(1000, num_spectra)
for chunk_start in 1:chunk_size:num_spectra
chunk_end = min(chunk_start + chunk_size - 1, num_spectra)
chunk_range = chunk_start:chunk_end
println("Processing chunk $chunk_start - $chunk_end / $num_spectra")
# Process current chunk
_iterate_spectra_fast(msi_data, collect(chunk_range)) do idx, mz, intensity
# Store metadata
modes[idx] = msi_data.spectra_metadata[idx].mode
if isempty(mz)
tics[idx] = 0.0
bpis[idx] = 0.0
bp_mzs[idx] = 0.0
num_points[idx] = 0
min_mzs[idx] = Inf
max_mzs[idx] = -Inf
return
end
# Update thread-local global m/z range
local_min, local_max = extrema(mz)
thread_id = Threads.threadid()
thread_local_min_mz[thread_id] = min(thread_local_min_mz[thread_id], local_min)
thread_local_max_mz[thread_id] = max(thread_local_max_mz[thread_id], local_max)
min_mzs[idx] = local_min
max_mzs[idx] = local_max
# Calculate per-spectrum stats
tics[idx] = sum(intensity)
max_int, max_idx = findmax(intensity)
bpis[idx] = max_int
bp_mzs[idx] = mz[max_idx]
num_points[idx] = length(mz)
end
# Force GC between chunks to control memory
if chunk_end % 5000 == 0
GC.gc()
end
end
# Combine thread-local global stats
g_min_mz = minimum(thread_local_min_mz)
g_max_mz = maximum(thread_local_max_mz)
# Create results
computed_stats_df = DataFrame(
SpectrumID = 1:num_spectra,
TIC = tics,
BPI = bpis,
BasePeakMZ = bp_mzs,
NumPoints = num_points,
MinMZ = min_mzs,
MaxMZ = max_mzs,
Mode = modes
)
# Set atomic values
Threads.atomic_xchg!(msi_data.global_min_mz, g_min_mz)
Threads.atomic_xchg!(msi_data.global_max_mz, g_max_mz)
# Set other data
msi_data.spectrum_stats_df = computed_stats_df
msi_data.bloom_filters = bloom_filters # Will be created on-demand
set!(msi_data.analytics_ready)
duration = (time_ns() - start_time) / 1e9
@printf "Streaming analytics complete in %.2f seconds.\n" duration
end
# Add mode statistics
centroid_count = count(==(CENTROID), modes)
profile_count = count(==(PROFILE), modes)
unknown_count = count(==(UNKNOWN), modes)
println("DEBUG Mode Statistics:")
println(" Centroid spectra: $centroid_count")
println(" Profile spectra: $profile_count")
println(" Unknown mode: $unknown_count")
# DEBUG: Print summary
println("DEBUG Analytics Summary:")
println(" Total spectra: $num_spectra")
println(" Spectra with data: $spectra_with_data")
println(" Empty spectra: $empty_spectra")
println(" Global m/z range: [$g_min_mz, $g_max_mz]")
println(" Centroid spectra: $(count(==(CENTROID), modes))")
println(" Profile spectra: $(count(==(PROFILE), modes))")
println(" Compressed spectra: $(sum(is_compressed))")
println(" Average points per spectrum: $(mean(num_points))")
# Populate the MSIData object
msi_data.global_min_mz = g_min_mz
msi_data.global_max_mz = g_max_mz
msi_data.spectrum_stats_df = DataFrame(
SpectrumID = 1:num_spectra,
TIC = tics,
BPI = bpis,
BasePeakMZ = bp_mzs,
NumPoints = num_points,
MinMZ = min_mzs,
MaxMZ = max_mzs,
Mode = modes,
IsCompressed = is_compressed
)
duration = (time_ns() - start_time) / 1e9
@printf "Analytics pre-computation complete in %.2f seconds.\n" duration
return
end
"""
get_total_spectrum_imzml(msi_data::MSIData; num_bins::Int=2000)
@ -661,10 +776,10 @@ function get_total_spectrum_imzml(msi_data::MSIData; num_bins::Int=2000, masked_
local global_min_mz, global_max_mz
if msi_data.global_min_mz !== nothing
if Threads.atomic_add!(msi_data.global_min_mz, 0.0) !== Inf
println(" Using pre-computed m/z range.")
global_min_mz = msi_data.global_min_mz
global_max_mz = msi_data.global_max_mz
global_min_mz = Threads.atomic_add!(msi_data.global_min_mz, 0.0)
global_max_mz = 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
pass1_start_time = time_ns()
@ -687,24 +802,26 @@ function get_total_spectrum_imzml(msi_data::MSIData; num_bins::Int=2000, masked_
# Use and cache the result
global_min_mz = g_min_mz
global_max_mz = g_max_mz
msi_data.global_min_mz = global_min_mz
msi_data.global_max_mz = global_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
# 2. Define Bins and precompute constants
mz_bins = range(global_min_mz, stop=global_max_mz, length=num_bins)
intensity_sum = zeros(Float64, num_bins)
bin_step = step(mz_bins)
inv_bin_step = 1.0 / bin_step # Precompute reciprocal to avoid division
min_mz = global_min_mz
# Initialize thread-local intensity sums
thread_local_intensity_sums = [zeros(Float64, num_bins) for _ in 1:Threads.nthreads()]
thread_local_spectra_processed = [0 for _ in 1:Threads.nthreads()]
# 3. Second Pass: Optimized binning
pass2_start_time = time_ns()
println(" Pass 2: Summing intensities into $num_bins bins...")
num_spectra_processed = 0
_iterate_spectra_fast(msi_data, masked_indices === nothing ? nothing : collect(masked_indices)) do idx, mz, intensity
num_spectra_processed += 1
thread_id = Threads.threadid()
thread_local_spectra_processed[thread_id] += 1
if isempty(mz)
return
end
@ -721,9 +838,17 @@ function get_total_spectrum_imzml(msi_data::MSIData; num_bins::Int=2000, masked_
# Branchless version using clamp
final_index = clamp(bin_index, 1, num_bins)
intensity_sum[final_index] += intensity_view[i]
thread_local_intensity_sums[thread_id][final_index] += intensity_view[i]
end
end
# Combine thread-local results
intensity_sum = zeros(Float64, num_bins)
num_spectra_processed = 0
for i in 1: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
@ -753,10 +878,10 @@ function get_total_spectrum_mzml(msi_data::MSIData; num_bins::Int=2000, masked_i
local global_min_mz, global_max_mz
if msi_data.global_min_mz !== nothing
if Threads.atomic_add!(msi_data.global_min_mz, 0.0) !== Inf
println(" Using pre-computed m/z range.")
global_min_mz = msi_data.global_min_mz
global_max_mz = msi_data.global_max_mz
global_min_mz = Threads.atomic_add!(msi_data.global_min_mz, 0.0)
global_max_mz = 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()
@ -780,8 +905,7 @@ function get_total_spectrum_mzml(msi_data::MSIData; num_bins::Int=2000, masked_i
# Use and cache the result
global_min_mz = g_min_mz
global_max_mz = g_max_mz
msi_data.global_min_mz = global_min_mz
msi_data.global_max_mz = global_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
@ -792,14 +916,17 @@ function get_total_spectrum_mzml(msi_data::MSIData; num_bins::Int=2000, masked_i
println(" Pass 2: Summing intensities into $num_bins bins...")
mz_bins = range(global_min_mz, stop=global_max_mz, length=num_bins)
intensity_sum = zeros(Float64, num_bins)
bin_step = step(mz_bins)
inv_bin_step = 1.0 / bin_step # Precompute reciprocal
min_mz = global_min_mz
num_spectra_processed = 0
# Initialize thread-local intensity sums
thread_local_intensity_sums = [zeros(Float64, num_bins) for _ in 1:Threads.nthreads()]
thread_local_spectra_processed = [0 for _ in 1:Threads.nthreads()]
_iterate_spectra_fast(msi_data, masked_indices === nothing ? nothing : collect(masked_indices)) do idx, mz, intensity
num_spectra_processed += 1
thread_id = Threads.threadid()
thread_local_spectra_processed[thread_id] += 1
if isempty(mz)
return
end
@ -811,9 +938,17 @@ function get_total_spectrum_mzml(msi_data::MSIData; num_bins::Int=2000, masked_i
# Branchless version using clamp
final_index = clamp(bin_index, 1, num_bins)
intensity_sum[final_index] += intensity[i]
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: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
@ -956,7 +1091,7 @@ the file stream `io` is already positioned at the correct offset.
- An error if zlib decompression fails, which can indicate corrupt data or
an incorrect offset in the `.imzML` metadata.
"""
function read_compressed_array(io::IO, asset::SpectrumAsset, format::Type)
function read_compressed_array(data::MSIData, io::IO, asset::SpectrumAsset, format::Type)
# Add validation before seeking
if asset.offset < 0 || asset.offset >= filesize(io)
throw(FileFormatError("Invalid asset offset: $(asset.offset) for file size $(filesize(io))"))
@ -964,30 +1099,43 @@ function read_compressed_array(io::IO, asset::SpectrumAsset, format::Type)
seek(io, asset.offset)
if asset.is_compressed
# Read compressed bytes - use encoded_length as the number of compressed bytes
compressed_bytes = read(io, asset.encoded_length)
# Get buffer for compressed bytes
compressed_bytes_buffer = get_buffer!(data.buffer_pool, asset.encoded_length)
readbytes!(io, compressed_bytes_buffer, asset.encoded_length)
println("DEBUG: Decompressing data - offset=$(asset.offset), compressed_bytes=$(length(compressed_bytes))")
println("DEBUG: Decompressing data - offset=$(asset.offset), compressed_bytes=$(length(compressed_bytes_buffer))")
local decompressed_bytes
local decompressed_bytes_buffer
try
decompressed_bytes = Libz.inflate(compressed_bytes)
println("DEBUG: Decompression successful - decompressed_bytes=$(length(decompressed_bytes))")
# Estimate decompressed size (can be larger than compressed)
# A common heuristic is 4x compressed size, but zlib can be more efficient
# For now, let Libz.inflate handle allocation, then copy to pooled buffer
# This is a temporary allocation, will be optimized later if needed
temp_decompressed = Libz.inflate(compressed_bytes_buffer)
decompressed_bytes_buffer = get_buffer!(data.buffer_pool, length(temp_decompressed))
copyto!(decompressed_bytes_buffer, temp_decompressed)
println("DEBUG: Decompression successful - decompressed_bytes=$(length(decompressed_bytes_buffer))")
catch e
@error "ZLIB DECOMPRESSION FAILED. This is likely due to an incorrect offset or corrupt data in the .ibd file."
@error "Asset offset: $(asset.offset), Encoded length: $(asset.encoded_length)"
# Print first 16 bytes to stderr for diagnosis
bytes_to_print = min(16, length(compressed_bytes))
bytes_to_print = min(16, length(compressed_bytes_buffer))
@error "First $bytes_to_print bytes of the data chunk we tried to decompress:"
println(stderr, view(compressed_bytes, 1:bytes_to_print))
println(stderr, view(compressed_bytes_buffer, 1:bytes_to_print))
rethrow(e)
finally
release_buffer!(data.buffer_pool, compressed_bytes_buffer)
end
# Use an IOBuffer to safely read the data
bytes_io = IOBuffer(decompressed_bytes)
bytes_io = IOBuffer(decompressed_bytes_buffer)
n_elements = bytes_io.size ÷ sizeof(format)
array = Array{format}(undef, n_elements)
read!(bytes_io, array)
release_buffer!(data.buffer_pool, decompressed_bytes_buffer)
return array
else
# Read uncompressed data directly
@ -1082,8 +1230,8 @@ function _iterate_compressed_fast(f::Function, data::MSIData, source::ImzMLSourc
end
# Read and decompress each array
mz_array = read_compressed_array(source.ibd_handle, meta.mz_asset, source.mz_format)
intensity_array = read_compressed_array(source.ibd_handle, meta.int_asset, source.intensity_format)
mz_array = read_compressed_array(data, source.ibd_handle, meta.mz_asset, source.mz_format)
intensity_array = read_compressed_array(data, source.ibd_handle, meta.int_asset, source.intensity_format)
mz_array .= ltoh.(mz_array)
intensity_array .= ltoh.(intensity_array)
@ -1143,15 +1291,15 @@ function _iterate_spectra_fast_impl(f::Function, data::MSIData, source::MzMLSour
for (i, _) in indices_with_offsets
meta = data.spectra_metadata[i]
mz = read_binary_vector(source.file_handle, meta.mz_asset)
intensity = read_binary_vector(source.file_handle, meta.int_asset)
mz = read_binary_vector(data, source.file_handle, meta.mz_asset)
intensity = read_binary_vector(data, source.file_handle, meta.int_asset)
f(i, mz, intensity)
end
end
"""
_iterate_spectra_fast(f::Function, data::MSIData, indices_to_iterate=nothing)
_iterate_spectra_fast_serial(f::Function, data::MSIData, indices_to_iterate=nothing)
An internal, high-performance iterator for bulk processing that bypasses the `GetSpectrum` cache.
It provides direct, sequential access to spectral data, making it the most efficient
@ -1173,7 +1321,183 @@ This is a low-level function that bypasses the public `GetSpectrum` API. It does
use the cache and is not guaranteed to be thread-safe if the same `MSIData` object
is accessed from multiple threads concurrently.
"""
function _iterate_spectra_fast(f::Function, data::MSIData, indices_to_iterate::Union{AbstractVector{Int}, Nothing}=nothing)
function _iterate_spectra_fast_serial(f::Function, data::MSIData, indices_to_iterate::Union{AbstractVector{Int}, Nothing}=nothing)
# Dispatch to the correct implementation based on the source type
_iterate_spectra_fast_impl(f, data, data.source, indices_to_iterate)
end
"""
_iterate_spectra_fast(f::Function, data::MSIData, indices_to_iterate=nothing)
A smart dispatcher for internal, high-performance iteration over spectra.
It automatically chooses between serial and parallel execution based on the
number of available threads (`Threads.nthreads()`).
# Arguments
- `f`: A function to call for each spectrum, with the signature `f(index, mz, intensity)`.
- `data`: The `MSIData` object.
- `indices_to_iterate`: An optional `AbstractVector{Int}` specifying a subset of
spectrum indices to iterate over. If `nothing`, it iterates over all spectra.
"""
function _iterate_spectra_fast(f::Function, data::MSIData, indices_to_iterate::Union{AbstractVector{Int}, Nothing}=nothing)
# Default to all indices if not specified
all_indices = (indices_to_iterate === nothing) ? (1:length(data.spectra_metadata)) : indices_to_iterate
if Threads.nthreads() > 1 && length(all_indices) > 100 # Heuristic: only parallelize for enough work
_iterate_spectra_fast_parallel(f, data, all_indices)
else
_iterate_spectra_fast_serial(f, data, all_indices)
end
end
"""
_iterate_spectra_fast_threadsafe(f::Function, data::MSIData, indices::AbstractVector{Int})
Thread-safe version of the fast iterator that uses the cache lock to ensure exclusive file access.
This is safe to use when multiple threads might access the same MSIData object concurrently.
# Arguments
- `f`: Function to call for each spectrum with signature `f(index, mz, intensity)`
- `data`: The MSIData object
- `indices`: Specific spectrum indices to iterate over
# Performance
- Slower than non-thread-safe version due to locking
- Safe for concurrent access to same MSIData object
- Good for small to medium datasets where parallel overhead isn't justified
"""
function _iterate_spectra_fast_threadsafe(f::Function, data::MSIData, indices::AbstractVector{Int})
# Use the cache lock to ensure exclusive access to file
lock(data.cache_lock) do
_iterate_spectra_fast(f, data, indices)
end
end
"""
_iterate_spectra_fast_parallel(f::Function, data::MSIData, indices::AbstractVector{Int})
Parallel version that creates thread-local file handles for maximum performance.
Each thread gets its own file handle, eliminating contention.
# Arguments
- `f`: Function to call for each spectrum with signature `f(index, mz, intensity)`
- `data`: The MSIData object
- `indices`: Specific spectrum indices to iterate over
# Performance
- Fastest for large datasets on multi-core systems
- Eliminates file I/O contention between threads
- Higher memory usage due to multiple file handles
- Best for bulk processing operations
"""
function _iterate_spectra_fast_parallel(f::Function, data::MSIData, indices::AbstractVector)
# Split indices into chunks for each thread
n_chunks = Threads.nthreads()
chunk_size = ceil(Int, length(indices) / n_chunks)
chunks = collect(Iterators.partition(indices, chunk_size))
Threads.@threads for chunk in chunks
# Each thread gets its own file handle based on source type
if data.source isa ImzMLSource
local_handle = open(data.source.ibd_handle.path, "r")
local_source = ImzMLSource(local_handle, data.source.mz_format, data.source.intensity_format)
try
# Use the appropriate implementation with thread-local source
_iterate_spectra_fast_impl(f, data, local_source, chunk)
finally
close(local_handle)
end
elseif data.source isa MzMLSource
local_handle = open(data.source.file_handle.path, "r")
local_source = MzMLSource(local_handle, data.source.mz_format, data.source.intensity_format)
try
_iterate_spectra_fast_impl(f, data, local_source, chunk)
finally
close(local_handle)
end
end
end
end
"""
_iterate_spectra_fast_parallel(f::Function, data::MSIData)
Parallel version that processes all spectra in the dataset.
Convenience wrapper for the indexed version.
"""
function _iterate_spectra_fast_parallel(f::Function, data::MSIData)
all_indices = collect(1:length(data.spectra_metadata))
_iterate_spectra_fast_parallel(f, data, all_indices)
end
"""
_iterate_spectra_fast_threadsafe(f::Function, data::MSIData)
Thread-safe version that processes all spectra in the dataset.
Convenience wrapper for the indexed version.
"""
function _iterate_spectra_fast_threadsafe(f::Function, data::MSIData)
all_indices = collect(1:length(data.spectra_metadata))
_iterate_spectra_fast_threadsafe(f, data, all_indices)
end
"""
get_bloom_filter(data::MSIData, index::Int)::BloomFilter{Int}
Retrieves the Bloom filter for a specific spectrum by its index.
If the Bloom filter has not yet been created (due to lazy initialization),
it will be created on-demand and cached.
# Arguments
- `data`: The `MSIData` object.
- `index`: The linear index of the spectrum.
# Returns
- A `BloomFilter{Int}` for the specified spectrum.
"""
function get_bloom_filter(data::MSIData, index::Int)::BloomFilter{Int}
lock(data.cache_lock) do
if data.bloom_filters === nothing
# Initialize Bloom filters on-demand
data.bloom_filters = Vector{Union{BloomFilter{Int}, Nothing}}(undef, length(data.spectra_metadata))
fill!(data.bloom_filters, nothing) # Initialize with nothing
end
if data.bloom_filters[index] === nothing
# Create Bloom filter lazily
mz, intensity = GetSpectrum(data, index)
data.bloom_filters[index] = create_bloom_filter_for_spectrum(mz)
end
return data.bloom_filters[index]
end
end
"""
monitor_memory_usage(phase::String)
Utility function to print current memory usage (allocated bytes and RSS)
at different phases of execution.
# Arguments
- `phase`: A string describing the current phase (e.g., "Initialization", "Processing chunk 1").
"""
function monitor_memory_usage(phase::String)
# Use Julia's internal memory tracking
bytes_allocated = Base.gc_bytes()
mb_allocated = bytes_allocated / 1024^2
# Get process memory from system
if Sys.islinux()
proc_status = read("/proc/self/status", String)
vm_rss_match = match(r"VmRSS:\s*(\d+)\s*kB", proc_status)
if vm_rss_match !== nothing
vm_rss_mb = parse(Int, vm_rss_match[1]) / 1024
@printf "[%s] Memory: %.2f MB allocated, %.2f MB RSS\n" phase mb_allocated vm_rss_mb
end
else
@printf "[%s] Memory: %.2f MB allocated\n" phase mb_allocated
end
end

6
src/MSI_src.jl
View File

@ -7,17 +7,16 @@ export OpenMSIData,
GetSpectrum,
IterateSpectra,
ImportMzmlFile,
load_slices,
plot_slices,
plot_slice,
get_total_spectrum,
get_average_spectrum,
LoadMzml,
precompute_analytics,
process_spectrum,
generate_colorbar_image,
process_image_pipeline,
load_and_prepare_mask
load_and_prepare_mask,
set_global_mz_range!
# Export the public Preprocessing API
export FeatureMatrix,
@ -40,6 +39,7 @@ export FeatureMatrix,
get_common_calibration_standards
# Include all source files directly into the main module
include("BloomFilters.jl")
include("Common.jl")
include("MSIData.jl")
include("ParserHelpers.jl")

491
src/imzML.jl
View File

@ -105,45 +105,6 @@ end
Reads metadata to determine the byte offsets and data types for reading spectra.
"""
function get_spectrum_attributes(stream::IO, hIbd::IO)
skip = Vector{UInt32}(undef, 8)
offset = get_spectrum_tag_offset(stream)
tag = find_tag(stream, r" accession=\"IMS:100005(\d)\"(.+)")
skip[1] = tag.captures[1][1] - '0' + 1
skip[2] = xor(skip[1], 3)
value = match(r" value=\"(\d+)\".+", tag.captures[2])
skip[5] = position(stream) - offset - length(value.match) - 2
value = match(r" value=\"\d+\".+", readline(stream))
skip[6] = value.offset
offset = position(stream)
tag = find_tag(stream, r"^\s*<referenceableParamGroupRef(.+)")
value = get_attribute(tag.captures[1], "ref")
skip[3] = (value.captures[1] == "intensityArray") + 3
skip[4] = xor(skip[3], 7)
k = 2
while k != 0
tag = find_tag(stream, r" accession=\"IMS:100010(\d)\"(.+)")
accession_val = tag.captures[1][1]
if accession_val == '2'
value = get_attribute(tag.match, "value")
seek(hIbd, parse(Int64, value.captures[1]))
k -= 1
elseif accession_val == '3'
skip[7] = position(stream) - offset - length(tag.match)
offset = position(stream)
k -= 1
end
end
find_tag(stream, r"^\s*</spectrum>")
skip[8] = position(stream) - offset
return skip
end
function get_spectrum_attributes_v2(stream::IO, hIbd::IO)
# Look for position x
readuntil(stream, "IMS:1000050")
x_skip = 0
@ -239,7 +200,7 @@ function load_imzml_lazy(file_path::String; cache_size::Int=100)
println("DEBUG: Opening file streams for .imzML and .ibd")
stream = open(file_path, "r")
hIbd = open(ibd_path, "r")
ts_hIbd = ThreadSafeFileHandle(ibd_path)
try
println("DEBUG: Configuring axes...")
@ -287,7 +248,7 @@ function load_imzml_lazy(file_path::String; cache_size::Int=100)
local spectra_metadata
if parser_type == :neofx
println("DEBUG: Using neofx parser.")
spectra_metadata = parse_neofx(stream, hIbd, param_groups, width, height, num_spectra,
spectra_metadata = parse_neofx(stream, ts_hIbd, param_groups, width, height, num_spectra,
default_mz_format, default_intensity_format,
mz_is_compressed, int_is_compressed, global_mode)
#=
@ -305,7 +266,7 @@ function load_imzml_lazy(file_path::String; cache_size::Int=100)
=#
else
println("DEBUG: Using compressed parser.")
spectra_metadata = parse_compressed(stream, hIbd, param_groups, width, height, num_spectra,
spectra_metadata = parse_compressed(stream, ts_hIbd, param_groups, width, height, num_spectra,
default_mz_format, default_intensity_format,
mz_is_compressed, int_is_compressed, global_mode)
end
@ -325,19 +286,23 @@ function load_imzml_lazy(file_path::String; cache_size::Int=100)
end
println("DEBUG: Coordinate map built.")
source = ImzMLSource(ts_hIbd, default_mz_format, default_intensity_format)
println("DEBUG: Creating MSIData object.")
msi_data = MSIData(source, spectra_metadata, (width, height), coordinate_map, cache_size)
# Close the XML stream as it's no longer needed
close(stream)
source = ImzMLSource(hIbd, default_mz_format, default_intensity_format)
println("DEBUG: Creating MSIData object.")
return MSIData(source, spectra_metadata, (width, height), coordinate_map, cache_size)
return msi_data
catch e
close(stream)
close(hIbd)
close(ts_hIbd) # Ensure IBD handle is closed on error
rethrow(e)
end
end
#=
function parse_uncompressed(stream::IO, hIbd::IO, param_groups::Dict{String, SpecDim},
width::Int32, height::Int32, num_spectra::Int32,
mz_format::DataType, intensity_format::DataType,
@ -350,7 +315,7 @@ function parse_uncompressed(stream::IO, hIbd::IO, param_groups::Dict{String, Spe
seek(hIbd, 16)
current_ibd_offset = position(hIbd)
attr = get_spectrum_attributes_v2(stream, hIbd)
attr = get_spectrum_attributes(stream, hIbd)
seek(stream, start_of_spectra_xml)
println("DEBUG: Initial IBD offset (after UUID): $current_ibd_offset")
@ -363,7 +328,7 @@ function parse_uncompressed(stream::IO, hIbd::IO, param_groups::Dict{String, Spe
# Find X coordinate
line = readuntil(stream, "IMS:1000050")
if eof(stream)
error("Unexpected EOF while looking for X coordinate")
throw(FileFormatError("Unexpected EOF while looking for X coordinate"))
end
val_tag_x = find_tag(stream, r"value=\"(\d+)\"")
x = parse(Int32, val_tag_x.captures[1])
@ -428,128 +393,139 @@ function parse_uncompressed(stream::IO, hIbd::IO, param_groups::Dict{String, Spe
return spectra_metadata
end
=#
function parse_compressed(stream::IO, hIbd::IO, param_groups::Dict{String, SpecDim},
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)
# Use concrete type for array_data
array_data_type = @NamedTuple{is_mz::Bool, array_length::Int32, encoded_length::Int64, offset::Int64}
for k in 1:num_spectra
# Read the full spectrum XML block
spectrum_buffer = IOBuffer()
# 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)
write(spectrum_buffer, line)
spectrum_start_line = line
break
end
end
# Parse lines within the spectrum block
while !eof(stream)
line = readline(stream)
write(spectrum_buffer, line)
if !occursin("<spectrum ", line) # Avoid re-reading the first line
line = readline(stream)
end
if occursin("</spectrum>", line)
break
end
end
spectrum_xml = String(take!(spectrum_buffer))
# Parse coordinates
x_match = match(r"IMS:1000050.*?value=\"(\d+)\"", spectrum_xml)
y_match = match(r"IMS:1000051.*?value=\"(\d+)\"", spectrum_xml)
x = x_match !== nothing ? parse(Int32, x_match.captures[1]) : Int32(0)
y = y_match !== nothing ? parse(Int32, y_match.captures[1]) : Int32(0)
# Parse mode - use SpectrumMode type directly
spectrum_mode = global_mode
if occursin("MS:1000127", spectrum_xml)
spectrum_mode = CENTROID
elseif occursin("MS:1000128", spectrum_xml)
spectrum_mode = PROFILE
end
# Parse binary data arrays with concrete typing
array_data = array_data_type[]
# Find all binaryDataArray blocks
array_matches = eachmatch(r"<binaryDataArray.*?<\/binaryDataArray>"s, spectrum_xml)
for array_match in array_matches
array_xml = array_match.match
# Determine if this is m/z or intensity array
is_mz = occursin("MS:1000514", array_xml) || occursin("mzArray", array_xml)
# Parse external data parameters
# Get array_length (nPoints)
array_len_cv_match = match(r"IMS:1000103.*?value=\"(\d+)\"", array_xml)
array_length = Int32(0)
if array_len_cv_match !== nothing
array_length = parse(Int32, array_len_cv_match.captures[1])
# Parse coordinates
x_match = match(r"IMS:1000050.*?value=\"(\d+)\"", line)
if x_match !== nothing
x = parse(Int32, x_match.captures[1])
end
if array_length == 0
nPoints_match = match(r"defaultArrayLength=\"(\d+)\"", spectrum_xml)
if nPoints_match !== nothing
array_length = parse(Int32, nPoints_match.captures[1])
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
# Get encoded_length
encoded_len_cv_match = match(r"IMS:1000104.*?value=\"(\d+)\"", array_xml)
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+)\"", array_xml)
if encoded_len_attr_match !== nothing
encoded_length = parse(Int64, encoded_len_attr_match.captures[1])
end
end
# Get offset
offset_match = match(r"IMS:1000102.*?value=\"(\d+)\"", array_xml)
offset = Int64(0)
if offset_match !== nothing
offset = parse(Int64, offset_match.captures[1])
end
if array_length > 0 && offset > 0
push!(array_data, (is_mz=is_mz, array_length=array_length,
encoded_length=encoded_length, offset=offset))
end
# Reset line to continue loop
line = ""
end
# Separate m/z and intensity arrays
mz_data = filter(d -> d.is_mz, array_data)
int_data = filter(d -> !d.is_mz, array_data)
# FIX: Handle empty spectra gracefully instead of throwing errors
current_ibd_pos = position(hIbd)
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
# Create placeholder metadata for empty/invalid spectrum
println("DEBUG: Spectrum $k is empty or invalid - creating placeholder metadata")
mz_asset = SpectrumAsset(default_mz_format, mz_is_compressed, current_ibd_pos, 0, :mz)
int_asset = SpectrumAsset(default_intensity_format, int_is_compressed, current_ibd_pos, 0, :intensity)
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]
# DEBUG: Print first spectrum details
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)")
println(" Expected m/z bytes: $(mz_info.array_length * sizeof(default_mz_format))")
println(" Expected intensity bytes: $(int_info.array_length * sizeof(default_intensity_format))")
end
# Create SpectrumAsset objects
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,
@ -562,7 +538,7 @@ function parse_compressed(stream::IO, hIbd::IO, param_groups::Dict{String, SpecD
return spectra_metadata
end
function parse_neofx(stream::IO, hIbd::IO, param_groups::Dict{String, SpecDim},
function parse_neofx(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)
@ -571,117 +547,128 @@ function parse_neofx(stream::IO, hIbd::IO, param_groups::Dict{String, SpecDim},
array_data_type = @NamedTuple{is_mz::Bool, array_length::Int32, encoded_length::Int64, offset::Int64}
for k in 1:num_spectra
# Read the full spectrum XML block
spectrum_buffer = IOBuffer()
# 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)
write(spectrum_buffer, line)
spectrum_start_line = line
break
end
end
# Parse lines within the spectrum block
while !eof(stream)
line = readline(stream)
write(spectrum_buffer, line)
if !occursin("<spectrum ", line) # Avoid re-reading the first line
line = readline(stream)
end
if occursin("</spectrum>", line)
break
end
end
spectrum_xml = String(take!(spectrum_buffer))
# Parse coordinates
x_match = match(r"IMS:1000050.*?value=\"(\d+)\"", spectrum_xml)
y_match = match(r"IMS:1000051.*?value=\"(\d+)\"", spectrum_xml)
x = x_match !== nothing ? parse(Int32, x_match.captures[1]) : Int32(0)
y = y_match !== nothing ? parse(Int32, y_match.captures[1]) : Int32(0)
# Parse mode
spectrum_mode = global_mode
if occursin("MS:1000127", spectrum_xml)
spectrum_mode = CENTROID
elseif occursin("MS:1000128", spectrum_xml)
spectrum_mode = PROFILE
end
# Parse binary data arrays
array_data = array_data_type[]
# Find all binaryDataArray blocks
array_matches = eachmatch(r"<binaryDataArray.*?<\/binaryDataArray>"s, spectrum_xml)
for array_match in array_matches
array_xml = array_match.match
# Determine if this is m/z or intensity array
is_mz = occursin("MS:1000514", array_xml) || occursin("mzArray", array_xml)
# Parse external data parameters
# Get array_length (nPoints)
array_len_cv_match = match(r"IMS:1000103.*?value=\"(\d+)\"", array_xml)
array_length = Int32(0)
if array_len_cv_match !== nothing
array_length = parse(Int32, array_len_cv_match.captures[1])
# Parse coordinates
x_match = match(r"IMS:1000050.*?value=\"(\d+)\"", line)
if x_match !== nothing
x = parse(Int32, x_match.captures[1])
end
if array_length == 0
nPoints_match = match(r"defaultArrayLength=\"(\d+)\"", spectrum_xml)
if nPoints_match !== nothing
array_length = parse(Int32, nPoints_match.captures[1])
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
# Get encoded_length
encoded_len_cv_match = match(r"IMS:1000104.*?value=\"(\d+)\"", array_xml)
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+)\"", array_xml)
if encoded_len_attr_match !== nothing
encoded_length = parse(Int64, encoded_len_attr_match.captures[1])
end
end
# Get offset
offset_match = match(r"IMS:1000102.*?value=\"(\d+)\"", array_xml)
offset = Int64(0)
if offset_match !== nothing
offset = parse(Int64, offset_match.captures[1])
end
if array_length > 0 && offset > 0
push!(array_data, (is_mz=is_mz, array_length=array_length,
encoded_length=encoded_length, offset=offset))
end
# Reset line to continue loop
line = ""
end
# Separate m/z and intensity arrays
mz_data = filter(d -> d.is_mz, array_data)
int_data = filter(d -> !d.is_mz, array_data)
# FIX: Handle empty spectra gracefully instead of throwing errors
current_ibd_pos = position(hIbd)
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
# Create placeholder metadata for empty/invalid spectrum
println("DEBUG: Spectrum $k is empty or invalid - creating placeholder metadata")
mz_asset = SpectrumAsset(default_mz_format, mz_is_compressed, current_ibd_pos, 0, :mz)
int_asset = SpectrumAsset(default_intensity_format, int_is_compressed, current_ibd_pos, 0, :intensity)
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]
# DEBUG: Print first spectrum details
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)")
println(" Expected m/z bytes: $(mz_info.array_length * sizeof(default_mz_format))")
println(" Expected intensity bytes: $(int_info.array_length * sizeof(default_intensity_format))")
end
# Create SpectrumAsset objects
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,
@ -739,6 +726,7 @@ function find_mass(mz_array::AbstractVector{<:Real}, intensity_array::AbstractVe
return max_intensity
end
#=
"""
load_slices(folder, masses, tolerance)
@ -795,7 +783,7 @@ function load_slices(folder::String, masses::AbstractVector{<:Real}, tolerance::
return (img_list, names)
end
=#
"""
get_mz_slice(data::MSIData, mass::Real, tolerance::Real)
@ -818,7 +806,7 @@ function get_mz_slice(data::MSIData, mass::Real, tolerance::Real; mask_path::Uni
end
# INTELLIGENT LOADING: Ensure analytics are computed for filtering.
if data.spectrum_stats_df === nothing || !hasproperty(data.spectrum_stats_df, :MinMZ)
if !is_set(data.analytics_ready)
println("Per-spectrum metadata not found. Running one-time analytics computation...")
precompute_analytics(data)
end
@ -826,13 +814,33 @@ function get_mz_slice(data::MSIData, mass::Real, tolerance::Real; mask_path::Uni
println("Using high-performance sequential iterator...")
target_min = mass - tolerance
target_max = mass + tolerance
stats_df = data.spectrum_stats_df
stats_df = get_spectrum_stats(data)
bloom_filters = get_bloom_filters(data)
# 1. Find all candidate spectra first for efficient filtering
candidate_indices = Set{Int}()
indices_to_check = masked_indices === nothing ? (1:length(data.spectra_metadata)) : masked_indices
for i in indices_to_check
# NEW: Bloom filter check with discretization
if bloom_filters !== nothing && !is_empty(bloom_filters[i])
discretization_factor = 100.0
min_mass_int = round(Int, (mass - tolerance) * discretization_factor)
max_mass_int = round(Int, (mass + tolerance) * discretization_factor)
found = false
for mass_int in min_mass_int:max_mass_int
if mass_int in bloom_filters[i]
found = true
break
end
end
if !found
continue # Definitely not in this spectrum
end
end
spec_min_mz = stats_df.MinMZ[i]
spec_max_mz = stats_df.MaxMZ[i]
if target_max >= spec_min_mz && target_min <= spec_max_mz
@ -890,13 +898,14 @@ function get_multiple_mz_slices(data::MSIData, masses::AbstractVector{<:Real}, t
end
# 2. Ensure analytics are computed for filtering.
if data.spectrum_stats_df === nothing || !hasproperty(data.spectrum_stats_df, :MinMZ)
if !is_set(data.analytics_ready)
println("Per-spectrum metadata not found. Running one-time analytics computation...")
precompute_analytics(data)
end
println("Filtering candidate spectra for $(length(masses)) m/z values...")
stats_df = data.spectrum_stats_df
stats_df = get_spectrum_stats(data)
bloom_filters = get_bloom_filters(data)
candidate_indices = Set{Int}()
indices_to_check = masked_indices === nothing ? (1:length(data.spectra_metadata)) : masked_indices
@ -909,6 +918,24 @@ function get_multiple_mz_slices(data::MSIData, masses::AbstractVector{<:Real}, t
if i in candidate_indices
continue
end
# NEW: Bloom filter check with discretization
if bloom_filters !== nothing && !is_empty(bloom_filters[i])
discretization_factor = 100.0
min_mass_int = round(Int, (mass - tolerance) * discretization_factor)
max_mass_int = round(Int, (mass + tolerance) * discretization_factor)
found = false
for mass_int in min_mass_int:max_mass_int
if mass_int in bloom_filters[i]
found = true
break
end
end
if !found
continue # Definitely not in this spectrum
end
end
spec_min_mz = stats_df.MinMZ[i]
spec_max_mz = stats_df.MaxMZ[i]
if target_max >= spec_min_mz && target_min <= spec_max_mz
@ -1272,7 +1299,7 @@ end
# Analysis and Visualization
#
# ============================================================================
#=
"""
display_statistics(slices::AbstractArray{<:AbstractMatrix{<:Real}, 2},
names::AbstractVector{String}, masses::AbstractVector{<:Real})
@ -1314,6 +1341,7 @@ function display_statistics(slices::AbstractArray{<:AbstractMatrix{<:Real}, 2},
end
return all_dfs
end
=#
"""
plot_slices(slices, names, masses, output_dir; stage_name, bins=256, dpi=150, global_bounds=nothing)
@ -1537,47 +1565,6 @@ function generate_palette(colorscheme, n_colors=256)
return palette
end
# ============================================================================
#
#
# Additional Performance Optimizations
#
# ============================================================================
"""
precompute_spectrum_bounds(msi_data::MSIData)
Precomputes min/max m/z values for each spectrum to enable fast filtering.
This avoids repeated computation during slice extraction.
"""
function precompute_spectrum_bounds(msi_data::MSIData)
n_spectra = length(msi_data.spectra_metadata)
min_mz = Vector{Float64}(undef, n_spectra)
max_mz = Vector{Float64}(undef, n_spectra)
_iterate_spectra_fast(msi_data) do idx, mz_array, intensity_array
min_mz[idx] = first(mz_array)
max_mz[idx] = last(mz_array)
end
return (min_mz, max_mz)
end
"""
process_columns_first!(A::AbstractMatrix)
Process a matrix column-by-column for better cache performance.
In Julia, arrays are stored in column-major order.
"""
function process_columns_first!(A::AbstractMatrix)
nrows, ncols = size(A)
@inbounds for col in 1:ncols, row in 1:nrows
# Process A[row, col] here
# This order is ~70% faster due to cache locality
end
return A
end
# Viridis color palette (256 colors) - defined as constant
const ViridisPalette = generate_palette(ColorSchemes.viridis)

40
src/mzML.jl
View File

@ -228,22 +228,22 @@ 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")
stream = open(file_path, "r")
ts_stream = ThreadSafeFileHandle(file_path, "r")
try
println("DEBUG: Finding index offset...")
index_offset = find_index_offset(stream)
index_offset = find_index_offset(ts_stream.handle)
println("DEBUG: Seeking to index list at offset $index_offset.")
seek(stream, index_offset)
seek(ts_stream.handle, index_offset)
println("DEBUG: Searching for '<index name=\"spectrum\">'.")
if find_tag(stream, 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(stream)
spectrum_offsets = parse_offset_list(ts_stream.handle)
if isempty(spectrum_offsets)
throw(FileFormatError("No spectrum offsets found."))
end
@ -256,7 +256,7 @@ function load_mzml_lazy(file_path::String; cache_size::Int=100)
# Use @inbounds for faster indexing in the loop
@inbounds for i in 1:num_spectra
spectra_metadata[i] = parse_spectrum_metadata(stream, spectrum_offsets[i])
spectra_metadata[i] = parse_spectrum_metadata(ts_stream.handle, spectrum_offsets[i])
# Progress reporting for large files
if i % 1000 == 0
@ -270,16 +270,17 @@ function load_mzml_lazy(file_path::String; cache_size::Int=100)
mz_format = first_meta.mz_asset.format
intensity_format = first_meta.int_asset.format
source = MzMLSource(stream, mz_format, intensity_format)
source = MzMLSource(ts_stream, mz_format, intensity_format)
println("DEBUG: Creating MSIData object.")
return MSIData(source, spectra_metadata, (0, 0), nothing, cache_size)
catch e
close(stream) # Ensure stream is closed on error
close(ts_stream) # Ensure stream is closed on error
rethrow(e)
end
end
#=
"""
LoadMzml(fileName::String)
@ -320,7 +321,9 @@ function LoadMzml(fileName::String)
return result_matrix
end
=#
#=
"""
load_mzml_batch(file_path::String, spectrum_indices::AbstractVector{Int})
@ -348,7 +351,9 @@ function load_mzml_batch(file_path::String, spectrum_indices::AbstractVector{Int
return results
end
=#
#=
"""
get_mzml_summary(file_path::String)::NamedTuple
@ -358,21 +363,21 @@ Quickly extracts summary information from an mzML file without loading all metad
- Named tuple with: num_spectra, mz_range, intensity_range, data_formats
"""
function get_mzml_summary(file_path::String)::NamedTuple
stream = open(file_path, "r")
ts_stream = ThreadSafeFileHandle(file_path, "r")
try
index_offset = find_index_offset(stream)
seek(stream, index_offset)
index_offset = find_index_offset(ts_stream.handle)
seek(ts_stream.handle, index_offset)
if find_tag(stream, r"<index\s+name=\"spectrum\"") === nothing
error("Could not find spectrum index.")
if find_tag(ts_stream.handle, r"<index\s+name=\"spectrum\"") === nothing
throw(FileFormatError("Could not find spectrum index."))
end
spectrum_offsets = parse_offset_list(stream)
spectrum_offsets = parse_offset_list(ts_stream.handle)
num_spectra = length(spectrum_offsets)
# Sample first few spectra to determine formats and ranges
sample_size = min(10, num_spectra)
sample_metadata = [parse_spectrum_metadata(stream, spectrum_offsets[i]) for i in 1:sample_size]
sample_metadata = [parse_spectrum_metadata(ts_stream.handle, spectrum_offsets[i]) for i in 1:sample_size]
# Extract formats from sample
mz_format = sample_metadata[1].mz_asset.format
@ -383,10 +388,12 @@ function get_mzml_summary(file_path::String)::NamedTuple
intensity_format=intensity_format,
sample_spectra=sample_size)
finally
close(stream)
close(ts_stream)
end
end
=#
#=
# Performance optimization: Cache frequently used regex patterns
const PRE_COMPILED_REGEX = (
encoded_length = r"<binaryDataArray\s+encodedLength=\"(\d+)\"",
@ -395,3 +402,4 @@ const PRE_COMPILED_REGEX = (
index_list = r"<index\s+name=\"spectrum\"",
index_offset = r"<indexListOffset>(\d+)</indexListOffset>"
)
=#

1
start_MSI_GUI.jl
View File

@ -4,6 +4,7 @@ Pkg.activate(".")
# Only instantiate in development mode
if get(ENV, "GENIE_ENV", "dev") != "prod"
@info "Development environment detected. Instantiating packages..."
Pkg.resolve()
Pkg.instantiate()
end
Pkg.gc()