diff --git a/Manifest.toml b/Manifest.toml index e1ae255..456263d 100644 --- a/Manifest.toml +++ b/Manifest.toml @@ -2,7 +2,7 @@ julia_version = "1.11.7" manifest_format = "2.0" -project_hash = "0f3959248d174e05d23b9e9267b00c6d2ead256d" +project_hash = "d4aecbadf6a54893101b3b89b57e8ae0ab391000" [[deps.ATK_jll]] deps = ["Artifacts", "Glib_jll", "JLLWrappers", "Libdl"] diff --git a/Project.toml b/Project.toml index 0db27f1..ff9fc31 100644 --- a/Project.toml +++ b/Project.toml @@ -15,6 +15,7 @@ GenieFramework = "a59fdf5c-6bf0-4f5d-949c-a137c9e2f353" Images = "916415d5-f1e6-5110-898d-aaa5f9f070e0" Libz = "2ec943e9-cfe8-584d-b93d-64dcb6d567b7" LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" +Mmap = "a63ad114-7e13-5084-954f-fe012c677804" NativeFileDialog = "e1fe445b-aa65-4df4-81c1-2041507f0fd4" NaturalSort = "c020b1a1-e9b0-503a-9c33-f039bfc54a85" PlotlyBase = "a03496cd-edff-5a9b-9e67-9cda94a718b5" diff --git a/app.jl b/app.jl index b836bce..f619549 100644 --- a/app.jl +++ b/app.jl @@ -15,6 +15,7 @@ using LinearAlgebra using NativeFileDialog # Opens the file explorer depending on the OS using StipplePlotly using Base.Filesystem: mv # To rename files in the system +using Printf # Required for @sprintf macro in colorbar generation # Bring MSIData into App module's scope using .MSI_src: MSIData, OpenMSIData, GetSpectrum, IterateSpectra, ImzMLSource, _iterate_spectra_fast, MzMLSource, find_mass, ViridisPalette, get_mz_slice, quantize_intensity, save_bitmap, median_filter, save_bitmap, downsample_spectrum, TrIQ @@ -295,16 +296,16 @@ include("./julia_imzML_visual.jl") fTime = time() eTime = round(fTime - sTime, digits=3) - msg = "File loaded in $(eTime) seconds. Calculating total spectrum..." + # msg = "File loaded in $(eTime) seconds. Calculating total spectrum..." + msg = "File loaded in $(eTime) seconds." # Enable UI controls btnStartDisable = !(msi_data.source isa ImzMLSource) btnPlotDisable = false btnSpectraDisable = false SpectraEnabled = true - + """ # --- Automatically generate and display the sum spectrum plot --- - # This is still part of the same async task progressSpectraPlot = true sTime = time() @@ -314,7 +315,7 @@ include("./julia_imzML_visual.jl") fTime = time() eTime = round(fTime - sTime, digits=3) msg = "Total spectrum plot loaded in $(eTime) seconds." - + """ catch e msi_data = nothing msg = "Error loading file: $e" @@ -324,7 +325,11 @@ include("./julia_imzML_visual.jl") SpectraEnabled = false @error "File loading failed" exception=(e, catch_backtrace()) finally - # This now correctly runs after everything is finished + # This block will always run at the end of the async task + GC.gc() + if Sys.islinux() + ccall(:malloc_trim, Int32, (Int32,), 0) + end progress = false progressSpectraPlot = false end @@ -347,7 +352,7 @@ include("./julia_imzML_visual.jl") msg = "No .imzML file loaded or selected file is not an .imzML. Please select a valid file." warning_msg = true elseif Nmass > 0 && Tol > 0 && Tol <= 1 && colorLevel > 1 && colorLevel < 257 - msg = "Creating image for Nmass=$(Nmass) Tol=$(Tol). Please be patient." + msg = "Creating image for m/z=$(Nmass) Tol=$(Tol). Please be patient." try # Use the new get_mz_slice with the centralized MSIData object slice = get_mz_slice(msi_data, Nmass, Tol) @@ -371,11 +376,8 @@ include("./julia_imzML_visual.jl") current_triq = "TrIQ_$(text_nmass).bmp" msgtriq = "TrIQ image with the Nmass of $(replace(text_nmass, "_" => "."))" - bound = MSI_src.get_outlier_thres(slice, triqProb) - levels = range(bound[1], stop=bound[2], length=8) - levels = vcat(levels, 2 * levels[end] - levels[end-1]) - Colorbar(fig[1, 1], colormap=cgrad(:viridis, colorLevel, categorical=true), limits=(0, bound[2]), ticks=levels, tickformat=log_tick_formatter, label="Intensity", size=25) - save("public/colorbar_TrIQ_$(text_nmass).png", fig) + colorbar_path = joinpath("public", "colorbar_TrIQ_$(text_nmass).png") + generate_colorbar_image(slice, colorLevel, colorbar_path, use_triq=true, triq_prob=triqProb) colorbarT = "/colorbar_TrIQ_$(text_nmass).png?t=$(timestamp)" current_col_triq = "colorbar_TrIQ_$(text_nmass).png" @@ -400,9 +402,8 @@ include("./julia_imzML_visual.jl") current_msi = "MSI_$(text_nmass).bmp" msgimg = "Image with the Nmass of $(replace(text_nmass, "_" => "."))" - levels = range(0, maximum(slice), length=8) - Colorbar(fig[1, 1], colormap=cgrad(:viridis, colorLevel, categorical=true), limits=(0, maximum(slice)), ticks=levels, tickformat=log_tick_formatter, label="Intensity", size=25) - save("public/colorbar_MSI_$(text_nmass).png", fig) + colorbar_path = joinpath("public", "colorbar_MSI_$(text_nmass).png") + generate_colorbar_image(slice, colorLevel, colorbar_path) colorbar = "/colorbar_MSI_$(text_nmass).png?t=$(timestamp)" current_col_msi = "colorbar_MSI_$(text_nmass).png" diff --git a/julia_imzML_visual.jl b/julia_imzML_visual.jl index e20abf7..666faff 100644 --- a/julia_imzML_visual.jl +++ b/julia_imzML_visual.jl @@ -325,7 +325,66 @@ function log_tick_formatter(values::Vector{Float64}) formValues[i]=value end return map((v, e) -> e == 0 ? "$(round(v, sigdigits=2))" : "$(round(v, sigdigits=2))x10" * Makie.UnicodeFun.to_superscript(e), formValues, exponents) +end +function generate_colorbar_image(slice_data::AbstractMatrix, color_levels::Int, output_path::String; use_triq::Bool=false, triq_prob::Float64=0.98) + # 1. Determine bounds based on whether TrIQ is used + min_val, max_val = if use_triq + MSI_src.get_outlier_thres(slice_data, triq_prob) + else + extrema(slice_data) + end + + # 2. Replicate the tick calculation logic from plot_slices + bins = color_levels + levels = range(min_val, stop=max_val, length=bins + 1) + level_range = levels[end] - levels[1] + + if level_range == 0 + levels = range(min_val - 0.1, stop=max_val + 0.1, length=bins + 1) + level_range = 0.2 + end + + exponent = level_range > 0 ? floor(log10(level_range)) / 3 : 0 + scale = 10^(3 * exponent) + scaled_levels = levels ./ scale + + format_num = level_range > 0 ? floor(log10(level_range)) % 3 : 0 + labels = if format_num == 0 + [ @sprintf("%3.2f", lvl) for lvl in scaled_levels] + elseif format_num == 1 + [ @sprintf("%3.2f", lvl) for lvl in scaled_levels] + else + [ @sprintf("%3.2f", lvl) for lvl in scaled_levels] + end + + divisors = 2:7 + remainders = (bins - 1) .% divisors + best_divisor = divisors[findlast(x -> x == minimum(remainders), remainders)] + tick_indices = round.(Int, range(1, stop=bins + 1, length=best_divisor + 1)) + + if !(1 in tick_indices) + pushfirst!(tick_indices, 1) + end + if !((bins + 1) in tick_indices) + push!(tick_indices, bins + 1) + end + unique!(sort!(tick_indices)) + + tick_positions = levels[tick_indices] + tick_labels = labels[tick_indices] + + # 3. Create and save the colorbar image + fig = Figure(size=(150, 250)) + Colorbar(fig[1, 1], + colormap=cgrad(:viridis, bins), + limits=(min_val, max_val), + label=(scale == 1 ? "Intensity" : "Intensity ×10^$(round(Int, 3 * exponent))"), + ticks=(tick_positions, tick_labels), + labelsize=20, + ticklabelsize=16 + ) + save(output_path, fig) end # meanSpectrumPlot recieves the local directory of the image as a string, @@ -408,7 +467,6 @@ function xySpectrumPlot(data::MSIData, xCoord::Int, yCoord::Int, imgWidth::Int, plotdata = [trace] plotlayout = layout - # Return the full data for other uses, and the plot data return plotdata, plotlayout, mz, intensity end diff --git a/src/MSIData.jl b/src/MSIData.jl index b570ee3..1bdfb61 100644 --- a/src/MSIData.jl +++ b/src/MSIData.jl @@ -6,27 +6,50 @@ including caching and iteration logic, for handling large mzML and imzML dataset efficiently. """ -using Base64, Libz, Serialization # For reading binary data +using Base64, Libz, Serialization, Printf # For reading binary data # Abstract type for different data sources (e.g., mzML, imzML) # This allows dispatching to the correct binary reading logic. abstract type MSDataSource end # Concrete type for imzML data sources +""" + ImzMLSource <: MSDataSource + +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 mz_format::Type intensity_format::Type end -# Concrete type for mzML data sources +""" + MzMLSource <: MSDataSource + +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 mz_format::Type intensity_format::Type end -# Struct to hold info about a binary data array (mz or intensity) +""" + SpectrumAsset + +Contains metadata for a single binary data array (m/z or intensity) within a spectrum. + +# Fields +- `format`: The data type of the elements (e.g., `Float32`, `Int64`). +- `is_compressed`: A boolean flag indicating if the data is compressed (e.g., with zlib). +- `offset`: The byte offset of the data within the file (`.ibd` for imzML, `.mzML` for mzML). +- `encoded_length`: The length of the data. For mzML, this is the Base64 encoded length. + For imzML, this is the number of elements in the array. +- `axis_type`: A symbol (`:mz` or `:intensity`) indicating the type of data. +""" struct SpectrumAsset format::Type is_compressed::Bool @@ -37,7 +60,17 @@ struct SpectrumAsset axis_type::Symbol end -# Metadata for a single spectrum, common to both formats +""" + SpectrumMetadata + +Contains all metadata for a single spectrum, common to both imzML and mzML formats. + +# Fields +- `x`, `y`: The spatial coordinates of the spectrum (for imzML only). +- `id`: The unique identifier string for the spectrum (for mzML only). +- `mz_asset`: A `SpectrumAsset` for the m/z array. +- `int_asset`: A `SpectrumAsset` for the intensity array. +""" struct SpectrumMetadata # For imzML x::Int32 @@ -51,7 +84,22 @@ struct SpectrumMetadata int_asset::SpectrumAsset end -# The main unified data object +""" + MSIData + +The primary object for interacting with mass spectrometry data. It provides a unified +interface for both `.mzML` and `.imzML` files and includes an LRU cache for +efficient repeated access to spectra. + +# Fields +- `source`: The underlying `MSDataSource` (`ImzMLSource` or `MzMLSource`). +- `spectra_metadata`: A vector of `SpectrumMetadata` for all spectra in the file. +- `image_dims`: A tuple `(width, height)` of the spatial dimensions (for imzML). +- `coordinate_map`: A matrix mapping `(x, y)` coordinates to a linear spectrum index (for imzML). +- `cache`: A dictionary holding cached spectra. +- `cache_order`: A vector tracking the usage order for the LRU cache. +- `cache_size`: The maximum number of spectra to store in the cache. +""" mutable struct MSIData source::MSDataSource spectra_metadata::Vector{SpectrumMetadata} @@ -81,6 +129,26 @@ end # --- Internal function for reading binary data --- # +""" + read_binary_vector(io::IO, asset::SpectrumAsset) + +Reads and decodes a single binary data vector (like m/z or intensity array) +from a `.mzML` file. The data is expected to be Base64-encoded and may be +compressed. + +This internal function handles: +1. Reading the raw Base64 string. +2. Decoding from Base64. +3. Decompressing the data if `asset.is_compressed` is true. +4. Converting the byte order from network (big-endian) to host order. + +# Arguments +- `io`: The IO stream of the `.mzML` file. +- `asset`: The `SpectrumAsset` containing metadata for the array. + +# Returns +- A `Vector` of the appropriate type containing the decoded data. +""" function read_binary_vector(io::IO, asset::SpectrumAsset) seek(io, asset.offset) raw_b64 = read(io, asset.encoded_length) @@ -94,6 +162,20 @@ function read_binary_vector(io::IO, asset::SpectrumAsset) end # Overload for different source types +""" + read_spectrum_from_disk(source::ImzMLSource, meta::SpectrumMetadata) + +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 +- `source`: An `ImzMLSource` containing the file handle and data formats. +- `meta`: The `SpectrumMetadata` for the spectrum to be read. + +# Returns +- A tuple `(mz, intensity)` containing the two requested data arrays. +""" function read_spectrum_from_disk(source::ImzMLSource, meta::SpectrumMetadata) # For imzML, the binary data is raw, not base64 encoded. # The `encoded_length` field in this case holds the number of points. @@ -106,9 +188,27 @@ function read_spectrum_from_disk(source::ImzMLSource, meta::SpectrumMetadata) seek(source.ibd_handle, meta.int_asset.offset) read!(source.ibd_handle, intensity) + # imzML data is little-endian. Convert to host byte order. + mz .= ltoh.(mz) + intensity .= ltoh.(intensity) + return mz, intensity end +""" + read_spectrum_from_disk(source::MzMLSource, meta::SpectrumMetadata) + +Reads a single spectrum's m/z and intensity arrays from a `.mzML` file. +This is a high-level function that orchestrates the reading process by calling +`read_binary_vector` for each of the m/z and intensity arrays. + +# Arguments +- `source`: A `MzMLSource` containing the file handle. +- `meta`: The `SpectrumMetadata` for the spectrum to be read. + +# Returns +- A tuple `(mz, intensity)` containing the two requested data arrays. +""" function read_spectrum_from_disk(source::MzMLSource, meta::SpectrumMetadata) # For mzML, data is Base64 encoded within the XML mz = read_binary_vector(source.file_handle, meta.mz_asset) @@ -182,13 +282,14 @@ end using Serialization function get_total_spectrum_imzml(msi_data::MSIData; num_bins::Int=2000) - println("Calculating total spectrum (2-pass method for imzML)...") + println("Calculating total spectrum for imzML (2-pass method)...") + total_start_time = time_ns() # 1. First Pass: Find the global m/z range + pass1_start_time = time_ns() println(" Pass 1: Finding global m/z range...") global_min_mz = Inf global_max_mz = -Inf - _iterate_spectra_fast(msi_data) do idx, mz, _ if !isempty(mz) local_min, local_max = extrema(mz) @@ -196,36 +297,66 @@ function get_total_spectrum_imzml(msi_data::MSIData; num_bins::Int=2000) global_max_mz = max(global_max_mz, local_max) end end + pass1_duration = (time_ns() - pass1_start_time) / 1e9 if !isfinite(global_min_mz) - @warn "Could not determine a valid m/z range. All spectra might be empty." + @warn "Could not determine a valid m/z range for imzML. All spectra might be empty." return (Float64[], Float64[]) end println(" Global m/z range found: [$(global_min_mz), $(global_max_mz)]") - # 2. Define Bins + # 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 - # 3. Second Pass: Sum intensities into bins + # 3. Second Pass: Optimized binning + pass2_start_time = time_ns() println(" Pass 2: Summing intensities into $num_bins bins...") _iterate_spectra_fast(msi_data) do idx, mz, intensity if isempty(mz) - return # continue equivalent + return end - for i in eachindex(mz) - bin_index = clamp(round(Int, (mz[i] - global_min_mz) / bin_step) + 1, 1, num_bins) - intensity_sum[bin_index] += intensity[i] + + # Use views to avoid bounds checks on the input arrays + mz_view = mz + intensity_view = intensity + + # Manual binning without clamp/round for SIMD + @inbounds for i in eachindex(mz_view) + # Calculate raw bin index (much faster than round+clamp) + raw_index = (mz_view[i] - min_mz) * inv_bin_step + 1.0 + bin_index = trunc(Int, raw_index) + + # Manual bounds checking (faster than clamp) + if 1 <= bin_index <= num_bins + intensity_sum[bin_index] += intensity_view[i] + elseif bin_index < 1 + intensity_sum[1] += intensity_view[i] + else # bin_index > num_bins + intensity_sum[num_bins] += intensity_view[i] + end end end - - println("Total spectrum calculation complete.") + pass2_duration = (time_ns() - pass2_start_time) / 1e9 + + total_duration = (time_ns() - total_start_time) / 1e9 + println("\n--- imzML Profiling ---") + @printf " Pass 1 (I/O only): %.2f seconds\n" pass1_duration + @printf " Pass 2 (I/O + Binning): %.2f seconds\n" pass2_duration + @printf " Est. Binning Overhead: %.2f seconds\n" (pass2_duration - pass1_duration) + @printf " Total Function Time: %.2f seconds\n" total_duration + println("-------------------------\n") + + println("Total spectrum calculation complete for imzML.") return (collect(mz_bins), intensity_sum) end function get_total_spectrum_mzml(msi_data::MSIData; num_bins::Int=2000) - println("Calculating total spectrum (single-pass optimization for mzML)...") + println("Calculating total spectrum for mzML (optimized single-pass method)...") + total_start_time = time_ns() num_spectra = length(msi_data.spectra_metadata) if num_spectra == 0 return (Float64[], Float64[]) @@ -233,31 +364,33 @@ function get_total_spectrum_mzml(msi_data::MSIData; num_bins::Int=2000) temp_path, temp_io = mktemp() try - # --- Pass 1: Write spectra to temp file and find min/max m/z --- - println(" Pass 1: Caching spectra and finding global m/z range...") + # --- Pass 1: Read from source, find m/z range, and write decoded spectra to temp file --- + pass1_start_time = time_ns() + println(" Pass 1: Caching decoded spectra and finding global m/z range...") global_min_mz = Inf global_max_mz = -Inf - + _iterate_spectra_fast(msi_data) do idx, mz, intensity - Serialization.serialize(temp_io, (mz, intensity)) if !isempty(mz) local_min, local_max = extrema(mz) global_min_mz = min(global_min_mz, local_min) global_max_mz = max(global_max_mz, local_max) end + Serialization.serialize(temp_io, (mz, intensity)) end flush(temp_io) + pass1_duration = (time_ns() - pass1_start_time) / 1e9 if !isfinite(global_min_mz) - @warn "Could not determine a valid m/z range. All spectra might be empty." + @warn "Could not determine a valid m/z range for mzML. All spectra might be empty." return (Float64[], Float64[]) end println(" Global m/z range found: [$(global_min_mz), $(global_max_mz)]") - # --- Pass 2: Read from temp file and bin intensities --- + # --- Pass 2: Read from fast temp file and bin intensities --- + pass2_start_time = time_ns() println(" Pass 2: Reading from cache and summing intensities into $num_bins bins...") - seekstart(temp_io) mz_bins = range(global_min_mz, stop=global_max_mz, length=num_bins) intensity_sum = zeros(Float64, num_bins) @@ -265,7 +398,6 @@ function get_total_spectrum_mzml(msi_data::MSIData; num_bins::Int=2000) while !eof(temp_io) mz, intensity = Serialization.deserialize(temp_io)::Tuple{AbstractVector, AbstractVector} - if isempty(mz) continue end @@ -274,18 +406,27 @@ function get_total_spectrum_mzml(msi_data::MSIData; num_bins::Int=2000) intensity_sum[bin_index] += intensity[i] end end - - println("Total spectrum calculation complete.") + pass2_duration = (time_ns() - pass2_start_time) / 1e9 + + total_duration = (time_ns() - total_start_time) / 1e9 + println("\n--- mzML Profiling ---") + @printf " Pass 1 (Read+Decode+Cache): %.2f seconds\n" pass1_duration + @printf " Pass 2 (Read Cache+Bin): %.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) finally close(temp_io) rm(temp_path, force=true) + println(" Temporary cache file removed.") end end """ - get_total_spectrum(msi_data::MSIData; num_bins::Int=20000) -> Tuple{Vector{Float64}, Vector{Float64}} + get_total_spectrum(msi_data::MSIData; num_bins::Int=2000) -> Tuple{Vector{Float64}, Vector{Float64}} Calculates the sum of all spectra in the dataset by binning. This function dispatches to a specialized implementation based on the file type @@ -301,7 +442,6 @@ function get_total_spectrum(msi_data::MSIData; num_bins::Int=2000) end end - """ get_average_spectrum(msi_data::MSIData; num_bins::Int=20000) -> Tuple{Vector{Float64}, Vector{Float64}} @@ -328,6 +468,12 @@ end # --- Iterator Implementation --- # +""" + MSIDataIterator + +A stateful iterator for sequentially accessing spectra from an `MSIData` object. +This is created by the `IterateSpectra` function. +""" struct MSIDataIterator data::MSIData end @@ -336,7 +482,7 @@ end IterateSpectra(data::MSIData) Returns an iterator that yields each spectrum, processing the file sequentially -with minimal memory overhead. +with minimal memory overhead. This iterator supports caching via `GetSpectrum`. This function is the core of the "Event-driven" access pattern. """ @@ -348,6 +494,20 @@ end Base.length(it::MSIDataIterator) = length(it.data.spectra_metadata) Base.eltype(it::MSIDataIterator) = Tuple{Int, Tuple{Vector, Vector}} # Index, (mz, intensity) +""" + Base.iterate(it::MSIDataIterator, state=1) + +Advances the iterator, returning the next spectrum in the sequence. It calls +`GetSpectrum`, which means the iteration process is cached according to the +`MSIData` object's cache settings. + +# Arguments +- `it`: The `MSIDataIterator` instance. +- `state`: The index of the next spectrum to retrieve. + +# Returns +- A tuple `((index, spectrum), next_state)` or `nothing` if iteration is complete. +""" function Base.iterate(it::MSIDataIterator, state=1) if state > length(it.data.spectra_metadata) return nothing # End of iteration @@ -364,35 +524,14 @@ end # --- High-performance Internal Iterator --- # """ - _iterate_spectra_fast(f::Function, data::MSIData) + _iterate_spectra_fast_impl(f::Function, data::MSIData, source::ImzMLSource) -An internal, high-performance iterator for bulk processing. -It avoids the overhead of `GetSpectrum` by reading directly from the file stream -and reusing pre-allocated buffers. - -- `f`: A function to call for each spectrum, with signature `f(index, mz, intensity)`. -- `data`: The `MSIData` object. +Internal implementation of the fast iterator for `.imzML` files. It reads +data directly from the `.ibd` file stream and reuses pre-allocated buffers +to minimize memory allocations and overhead, making it ideal for bulk processing tasks. """ function _iterate_spectra_fast_impl(f::Function, data::MSIData, source::ImzMLSource) - # This implementation is for imzML and is optimized for performance by - # using pre-allocated buffers. - max_points = 0 - for meta in data.spectra_metadata - # For imzML, encoded_length is the number of points - max_points = max(max_points, meta.mz_asset.encoded_length) - end - - # If all spectra are empty, just call f with empty arrays - if max_points == 0 && !isempty(data.spectra_metadata) - for i in 1:length(data.spectra_metadata) - f(i, source.mz_format[], source.intensity_format[]) - end - return - end - - mz_buffer = Vector{source.mz_format}(undef, max_points) - int_buffer = Vector{source.intensity_format}(undef, max_points) - + # Optimized implementation: allocates arrays per spectrum and uses an efficient I/O pattern. for i in 1:length(data.spectra_metadata) meta = data.spectra_metadata[i] nPoints = meta.mz_asset.encoded_length @@ -402,19 +541,43 @@ function _iterate_spectra_fast_impl(f::Function, data::MSIData, source::ImzMLSou continue end - mz_view = view(mz_buffer, 1:nPoints) - int_view = view(int_buffer, 1:nPoints) + # Allocate fresh arrays for each spectrum. + mz = Vector{source.mz_format}(undef, nPoints) + intensity = Vector{source.intensity_format}(undef, nPoints) - seek(source.ibd_handle, meta.mz_asset.offset) - read!(source.ibd_handle, mz_view) + # --- I/O OPTIMIZATION: Seek only once per spectrum --- + # The mz and intensity data are contiguous, so after reading the first, + # we can immediately read the second without a costly second seek. + if meta.mz_asset.offset < meta.int_asset.offset + # m/z is first, seek to it. + seek(source.ibd_handle, meta.mz_asset.offset) + # Read m/z, then immediately read intensity. + read!(source.ibd_handle, mz) + read!(source.ibd_handle, intensity) + else + # intensity is first, seek to it. + seek(source.ibd_handle, meta.int_asset.offset) + # Read intensity, then immediately read m/z. + read!(source.ibd_handle, intensity) + read!(source.ibd_handle, mz) + end - seek(source.ibd_handle, meta.int_asset.offset) - read!(source.ibd_handle, int_view) + # Convert byte order. + mz .= ltoh.(mz) + intensity .= ltoh.(intensity) - f(i, mz_view, int_view) + f(i, mz, intensity) end end +""" + _iterate_spectra_fast_impl(f::Function, data::MSIData, source::MzMLSource) + +Internal implementation of the fast iterator for `.mzML` files. It iterates +through each spectrum, decodes the Base64 data on the fly, and calls the +provided function. It bypasses the `GetSpectrum` cache to avoid storing +all decoded spectra in memory. +""" function _iterate_spectra_fast_impl(f::Function, data::MSIData, source::MzMLSource) # This implementation is for mzML. It iterates through each spectrum, # decodes it, and then calls the function `f`. It's less performant than @@ -430,6 +593,16 @@ function _iterate_spectra_fast_impl(f::Function, data::MSIData, source::MzMLSour end end +""" + _iterate_spectra_fast(f::Function, data::MSIData) + +An internal, high-performance iterator for bulk processing that bypasses the cache. +It dispatches to a specialized implementation based on the data source type +(`ImzMLSource` or `MzMLSource`). + +- `f`: A function to call for each spectrum, with signature `f(index, mz, intensity)`. +- `data`: The `MSIData` object. +""" function _iterate_spectra_fast(f::Function, data::MSIData) # Dispatch to the correct implementation based on the source type _iterate_spectra_fast_impl(f, data, data.source) diff --git a/src/MSI_src.jl b/src/MSI_src.jl index c6f3ce0..59a40c4 100644 --- a/src/MSI_src.jl +++ b/src/MSI_src.jl @@ -20,7 +20,7 @@ Opens a .mzML or .imzML file and prepares it for data access. This is the main entry point for the new data access API. """ -function OpenMSIData(filepath::String; cache_size=100) +function OpenMSIData(filepath::String; cache_size=300) if endswith(lowercase(filepath), ".mzml") return load_mzml_lazy(filepath, cache_size=cache_size) elseif endswith(lowercase(filepath), ".imzml") diff --git a/src/MzmlConverter.jl b/src/MzmlConverter.jl index 695eed4..f50ec52 100644 --- a/src/MzmlConverter.jl +++ b/src/MzmlConverter.jl @@ -11,13 +11,33 @@ using DataFrames, Printf, CSV # A struct to hold the processed pixel data before exporting +""" + ProcessedPixel + +A temporary struct to hold the data for a single, fully rendered pixel before +it is written to the final `.ibd` file. + +# Fields +- `coords`: A tuple `(x, y)` of the pixel's spatial coordinates. +- `mz`: The m/z array for the pixel. +- `intensity`: The calculated intensity array for the pixel. +""" struct ProcessedPixel coords::Tuple{Int, Int} mz::Vector{Float64} # Assuming m/z is consistent, can be optimized later intensity::Vector{Float32} end -# This struct will hold the metadata needed for the XML file +""" + BinaryMetadata + +A struct to hold the file offset and length for a spectrum's binary data arrays +(m/z and intensity) after they have been written to the `.ibd` file. + +# Fields +- `mz_offset`, `mz_length`: Byte offset and length for the m/z array. +- `int_offset`, `int_length`: Byte offset and length for the intensity array. +""" struct BinaryMetadata mz_offset::UInt64 mz_length::UInt64 @@ -28,8 +48,15 @@ end """ GetMzmlScanTime_linebyline(fileName::String) -(Internal) Slow, line-by-line parser for scan times. Used as a fallback if -the .mzML file index is missing or corrupt. +Parses a `.mzML` file line-by-line to extract the scan start time for each +spectrum. This is a slow and memory-intensive fallback method used only when +the faster, index-based `GetMzmlScanTime` fails. + +# Arguments +- `fileName`: Path to the `.mzML` file. + +# Returns +- A `Matrix{Int64}` where each row is `[spectrum_index, time_in_milliseconds]`. """ function GetMzmlScanTime_linebyline(fileName::String) times = Tuple{Int64, Int64}[] @@ -350,6 +377,28 @@ function MatchAcquireTime(sync_file_path::String, scans::Matrix{Int64}; img_widt end +""" + RenderPixel(pixel_info, scans, msi_data, scan_time_deltas, pixel_time_deltas) + +Reconstructs the spectrum for a single pixel by combining intensities from the +raw MS scans that occurred during the pixel's acquisition time. It uses a +weighted interpolation scheme based on the relative timing of scans and pixels. + +This function handles three cases: +1. A single scan falls entirely within the pixel's time window. +2. The pixel's time window is covered by two partial scans. +3. The pixel's time window covers one or more full scans plus two partial scans. + +# Arguments +- `pixel_info`: A row from the timing matrix containing the pixel's time and scan indices. +- `scans`: The matrix of scan times. +- `msi_data`: The `MSIData` object for the source `.mzML` file. +- `scan_time_deltas`: Pre-calculated time durations for each scan. +- `pixel_time_deltas`: Pre-calculated time durations for each pixel. + +# Returns +- A tuple `(mz_array, intensity_array)` for the rendered pixel spectrum. +""" function RenderPixel(pixel_info, scans, msi_data::MSIData, scan_time_deltas, pixel_time_deltas) pixel_time = pixel_info[3] first_scan = pixel_info[4] @@ -408,6 +457,26 @@ function RenderPixel(pixel_info, scans, msi_data::MSIData, scan_time_deltas, pix return (mz_array, new_intensity) end +""" + ConvertMzmlToImzml(source_file, target_ibd_file, timing_matrix, scans) + +Orchestrates the conversion of spectra from a `.mzML` file into a binary `.ibd` +file. It iterates through each pixel defined in the `timing_matrix`, calls +`RenderPixel` to reconstruct the pixel's spectrum, and writes the resulting +m/z and intensity arrays to the `.ibd` file in little-endian byte order. + +# Arguments +- `source_file`: Path to the source `.mzML` file. +- `target_ibd_file`: Path for the output `.ibd` binary file. +- `timing_matrix`: The output from `MatchAcquireTime`, mapping pixels to scans. +- `scans`: The matrix of scan times from `GetMzmlScanTime`. + +# Returns +- A tuple `(binary_meta_vec, coords_vec, (width, height))` containing: + - A vector of `BinaryMetadata` for each spectrum. + - A vector of `(x, y)` coordinate tuples. + - A tuple of the final image dimensions. +""" function ConvertMzmlToImzml(source_file::String, target_ibd_file::String, timing_matrix::Matrix{Int64}, scans::Matrix{Int64}) if size(timing_matrix, 1) == 0 # Create an empty .ibd file if there's nothing to process @@ -489,6 +558,23 @@ function ConvertMzmlToImzml(source_file::String, target_ibd_file::String, timing return binary_meta_vec, coords_vec, (width, height) end +""" + ExportImzml(target_file, binary_meta, coords, dims) + +Generates the `.imzML` metadata file. This XML file contains all the necessary +metadata to interpret the corresponding `.ibd` binary file, including references +to external data offsets, image dimensions, and CV parameters describing the +experiment and data format. + +# Arguments +- `target_file`: The path for the output `.imzML` file. +- `binary_meta`: A vector of `BinaryMetadata` structs with offset and length info. +- `coords`: A vector of `(x, y)` coordinates for each spectrum. +- `dims`: A tuple `(width, height)` of the final image dimensions. + +# Returns +- `true` on success, `false` on failure. +""" function ExportImzml(target_file::String, binary_meta::Vector{BinaryMetadata}, coords::Vector{Tuple{Int, Int}}, dims::Tuple{Int, Int}) ibd_file = replace(target_file, r"\.imzML$"i => ".ibd") diff --git a/src/ParserHelpers.jl b/src/ParserHelpers.jl index e0b21d6..5f52b95 100644 --- a/src/ParserHelpers.jl +++ b/src/ParserHelpers.jl @@ -59,7 +59,15 @@ end """ configure_spec_dim(stream) -Fills a `SpecDim` struct by parsing `cvParam` tags from the stream. +Parses a block of `` tags from an XML stream to configure a `SpecDim` +struct. It reads accessions to determine the data format (e.g., `Float32`), +compression status (`zlib`), and axis type (m/z vs. intensity). + +# Arguments +- `stream`: An IO stream positioned at the start of the `cvParam` block. + +# Returns +- A `SpecDim` struct populated with the parsed configuration. """ function configure_spec_dim(stream) axis = SpecDim(Float64, false, 1, 0) @@ -102,6 +110,18 @@ end # # ============================================================================ +""" + CVParams + +A mutable struct used during the parsing of mzML files to temporarily hold +CV parameter information for a data array before it is converted into a +`SpectrumAsset`. + +# Fields +- `format`: The data type (e.g., `Float64`). +- `is_compressed`: A boolean indicating if the data is compressed. +- `axis_type`: A symbol indicating the array type (`:mz` or `:intensity`). +""" mutable struct CVParams format::Type is_compressed::Bool diff --git a/src/imzML.jl b/src/imzML.jl index 0426109..61f60c3 100644 --- a/src/imzML.jl +++ b/src/imzML.jl @@ -1,3 +1,4 @@ +# src/imzML.jl using Images, Statistics, CairoMakie, DataFrames, Printf, ColorSchemes, StatsBase # --- Extracted from imzML.jl --- @@ -171,7 +172,10 @@ function load_imzml_lazy(file_path::String; cache_size=100) int_config_idx = findfirst(a -> a.Axis == 2, axis) mz_format = axis[mz_config_idx].Format intensity_format = axis[int_config_idx].Format + mz_is_compressed = axis[mz_config_idx].Packed + int_is_compressed = axis[int_config_idx].Packed println("DEBUG: m/z format: $mz_format, Intensity format: $intensity_format") + println("DEBUG: m/z compressed: $mz_is_compressed, Intensity compressed: $int_is_compressed") # --- NEW PARSING LOGIC based on the old, working code --- println("DEBUG: Learning file structure from first spectrum...") @@ -212,8 +216,8 @@ function load_imzml_lazy(file_path::String; cache_size=100) end # Create modern SpectrumAsset objects - mz_asset = SpectrumAsset(mz_format, false, mz_offset, nPoints, :mz) - int_asset = SpectrumAsset(intensity_format, false, int_offset, nPoints, :intensity) + mz_asset = SpectrumAsset(mz_format, mz_is_compressed, mz_offset, nPoints, :mz) + int_asset = SpectrumAsset(intensity_format, int_is_compressed, int_offset, nPoints, :intensity) spectra_metadata[k] = SpectrumMetadata(x, y, "", mz_asset, int_asset) @@ -449,6 +453,21 @@ end # # ============================================================================ +""" + get_outlier_thres(img, prob=0.98) + +Calculates dynamic intensity range bounds for an image based on a cumulative +probability histogram. This function replicates an R algorithm to find an +intensity threshold that corresponds to a given cumulative probability, which +is used to exclude outliers before normalization. + +# Arguments +- `img`: The input image matrix. +- `prob`: The cumulative probability threshold (default: 0.98) for outlier detection. + +# Returns +- A tuple `(low, high)` representing the calculated lower and upper intensity bounds. +""" function get_outlier_thres(img, prob=0.98) # DO NOT filter zeros. Use all pixel values like R does. int_values = vec(img) @@ -493,6 +512,21 @@ function get_outlier_thres(img, prob=0.98) return (low, actual_threshold) end +""" + set_pixel_depth(img, bounds, depth) + +Quantizes the intensity values of an image into a specified number of bins +(`depth`) within a given intensity range (`bounds`). Pixels outside the bounds +are clipped. + +# Arguments +- `img`: The input image matrix. +- `bounds`: A tuple `(min, max)` specifying the intensity range for quantization. +- `depth`: The number of quantization levels (bins). + +# Returns +- A `Matrix{UInt8}` with pixel values quantized to the specified depth. +""" function set_pixel_depth(img, bounds, depth) min_val, max_val = bounds bins = depth - 1 @@ -595,6 +629,19 @@ function quantize_intensity(slice::AbstractMatrix{<:Real}, levels::Integer=256) return image end +""" + median_filter(img) + +Applies a 3x3 median filter to the input image. This is a simple noise +reduction technique that replaces each pixel's value with the median value of +its 3x3 neighborhood. + +# Arguments +- `img`: The input image matrix. + +# Returns +- A new matrix containing the filtered image. +""" function median_filter(img) # 3x3 median filter implementation return mapwindow(median, img, (3, 3)) @@ -685,6 +732,28 @@ function display_statistics(slices, names, masses) return all_dfs end +""" + plot_slices(slices, names, masses, output_dir; stage_name, bins=256, dpi=150, global_bounds=nothing) + +Generates and saves a grid of image slice plots. Each row corresponds to a +file and each column to a mass, creating a comprehensive overview. + +# Arguments +- `slices`: A 2D array of image slice matrices (`n_files` x `n_masses`). +- `names`: A vector of file names, used for titling rows. +- `masses`: A vector of m/z values, used for titling columns. +- `output_dir`: The directory where the output plots will be saved. + +# Keyword Arguments +- `stage_name`: A string used to name the output files (e.g., "raw", "normalized"). +- `bins`: The number of color levels in the heatmap palette. +- `dpi`: The resolution for the saved image files. +- `global_bounds`: A vector of `(min, max)` tuples, one for each mass, to ensure a consistent + color scale across all files for a given mass. If `nothing`, bounds are calculated automatically. + +# Returns +- The generated `Figure` object from Makie. +""" function plot_slices(slices, names, masses, output_dir; stage_name, bins=256, dpi=150, global_bounds=nothing) n_files, n_masses = size(slices) diff --git a/src/mzML.jl b/src/mzML.jl index d65395a..90ad8f2 100644 --- a/src/mzML.jl +++ b/src/mzML.jl @@ -1,9 +1,22 @@ -# /home/pixel/Documents/Cinvestav_2025/JuliaMSI/src/mzML.jl +# src/mzML.jl # This file is responsible for parsing metadata from .mzML files. # It has been refactored to produce a unified MSIData object. -# This function is kept internal to the mzML parsing process +""" + get_spectrum_asset_metadata(stream) + +Parses a `` block within an mzML file to extract metadata +for a single data array (e.g., m/z or intensity). It reads CV parameters to +determine the data type, compression, and axis type. + +# Arguments +- `stream`: An IO stream positioned at the beginning of a `` block. + +# Returns +- A `SpectrumAsset` struct containing the parsed metadata, including the binary + data offset, encoded length, format, and compression status. +""" function get_spectrum_asset_metadata(stream) start_pos = position(stream) @@ -56,6 +69,10 @@ function get_spectrum_asset_metadata(stream) elseif acc_str == "MS:1000523" data_format = Float64 # println("DEBUG: Set format to Float64") + # MS:1000572 is a parent term for compression. While not ideal, if present, assume compression. + elseif acc_str == "MS:1000572" + compression_flag = true + # println("DEBUG: Set is_compressed to true based on parent term") elseif acc_str == "MS:1000574" compression_flag = true # println("DEBUG: Set is_compressed to true") @@ -80,6 +97,20 @@ function get_spectrum_asset_metadata(stream) end # This function is updated to return the generic SpectrumMetadata struct +""" + parse_spectrum_metadata(stream, offset::Int64) + +Parses an entire `` block from an mzML file, given a starting offset. +It extracts the spectrum ID and calls `get_spectrum_asset_metadata` to parse +the m/z and intensity array metadata. + +# Arguments +- `stream`: An IO stream for the mzML file. +- `offset`: The byte offset where the `` block begins. + +# Returns +- A `SpectrumMetadata` struct containing the parsed metadata for one spectrum. +""" function parse_spectrum_metadata(stream, offset::Int64) seek(stream, offset) @@ -97,6 +128,18 @@ function parse_spectrum_metadata(stream, offset::Int64) return SpectrumMetadata(0, 0, id, mz_asset, int_asset) end +""" + parse_offset_list(stream) + +Parses the `` block in an indexed mzML file to extract +the byte offsets for each spectrum. + +# Arguments +- `stream`: An IO stream positioned at the start of the `` block. + +# Returns +- A `Vector{Int64}` containing the byte offsets for all spectra. +""" function parse_offset_list(stream) offsets = Int64[] offset_regex = r"]*>(\d+)" @@ -121,6 +164,20 @@ function parse_offset_list(stream) end # This is the main lazy-loading function for mzML, now returning an MSIData object. +""" + load_mzml_lazy(file_path::String; cache_size::Int=100) + +Lazily loads an indexed `.mzML` file by parsing only the metadata. It reads the +spectrum index from the end of the file to get the offsets of each spectrum, +then parses the metadata for each spectrum without loading the binary data. + +# Arguments +- `file_path`: The path to the `.mzML` file. +- `cache_size`: The number of spectra to hold in an LRU cache for faster access. + +# Returns +- An `MSIData` object ready for lazy data access. +""" function load_mzml_lazy(file_path::String; cache_size::Int=100) println("DEBUG: Opening file stream for $file_path") stream = open(file_path, "r") @@ -181,8 +238,16 @@ end """ LoadMzml(fileName::String) -Eagerly loads all spectra from a .mzML file. -Provided for backward compatibility. Now uses the new MSIData architecture. +Eagerly loads all spectra from a .mzML file into memory. +This function now uses the new lazy-loading +architecture internally but presents the data in the old format. + +# Arguments +- `fileName`: The path to the `.mzML` file. + +# Returns +- A `2xN` matrix where `N` is the number of spectra. The first row contains + m/z arrays and the second row contains intensity arrays. """ function LoadMzml(fileName::String) # Use the lazy loader to get the MSIData object diff --git a/test/run_tests.jl b/test/run_tests.jl index 9b2efb2..d8c9de9 100644 --- a/test/run_tests.jl +++ b/test/run_tests.jl @@ -31,16 +31,18 @@ using MSI_src # --- Test Case 1: Standard .mzML file --- # A regular, non-imaging mzML file. # const TEST_MZML_FILE = "/home/pixel/Documents/Cinvestav_2025/Analisis/mzML/T9_A1.mzML" -# const TEST_MZML_FILE = "/home/pixel/Documents/Cinvestav_2025/Analisis/CE4_BF_R1/CE4_BF_R1.mzML" -# const TEST_MZML_FILE = "/home/pixel/Documents/Cinvestav_2025/Analisis/CE4_BF_R1/CE4_BF_R1.mzML" -const TEST_MZML_FILE = "/home/pixel/Documents/Cinvestav_2025/Analisis/Imaging_paper_spray/Imaging_paper_spray.mzML" +const TEST_MZML_FILE = "/home/pixel/Documents/Cinvestav_2025/Analisis/CE4_BF_R1/CE4_BF_R1.mzML" +# const TEST_MZML_FILE = "/home/pixel/Documents/Cinvestav_2025/Analisis/Imaging_paper_spray/Imaging_paper_spray.mzML" +# const TEST_MZML_FILE = "/home/pixel/Documents/Cinvestav_2025/Analisis/Imaging prueba Roya 1/Roya.mzML" const SPECTRUM_TO_PLOT = 1 # Which spectrum to plot from the file # --- Test Case 2: .mzML + Sync File for Conversion --- # The special .mzML file with one spectrum per pixel. -const CONVERSION_SOURCE_MZML = "/home/pixel/Documents/Cinvestav_2025/Analisis/CE4_BF_R1/CE4_BF_R1.mzML" +# const CONVERSION_SOURCE_MZML = "/home/pixel/Documents/Cinvestav_2025/Analisis/CE4_BF_R1/CE4_BF_R1.mzML" +const CONVERSION_SOURCE_MZML = TEST_MZML_FILE # The corresponding synchronization text file. const CONVERSION_SYNC_FILE = "/home/pixel/Documents/Cinvestav_2025/Analisis/CE4_BF_R1/CE4_BF_R1.txt" +# const CONVERSION_SYNC_FILE = "/home/pixel/Documents/Cinvestav_2025/Analisis/Imaging prueba Roya 1/Synchro.txt" # const CONVERSION_SOURCE_MZML = "/home/pixel/Documents/Cinvestav_2025/Analisis/Imaging_paper_spray/Imaging_paper_spray.mzML" # const CONVERSION_SYNC_FILE = "/home/pixel/Documents/Cinvestav_2025/Analisis/Imaging_paper_spray/Imaging_paper_spray.txt" @@ -50,18 +52,18 @@ const CONVERSION_TARGET_IMZML = "test/results/converted_mzml.imzML" # --- Test Case 3: Standard .imzML file --- # An existing imzML file (can be the one generated from Case 2). -# const TEST_IMZML_FILE = CONVERSION_TARGET_IMZML # The output from case 2 +const TEST_IMZML_FILE = CONVERSION_TARGET_IMZML # The output from case 2 # const TEST_IMZML_FILE = "/home/pixel/Documents/Cinvestav_2025/Analisis/imzML_AP_SMALDI/HR2MSImouseurinarybladderS096.imzML" # const TEST_IMZML_FILE = "/home/pixel/Documents/Cinvestav_2025/Analisis/Imaging_paper_spray/Imaging_paper_spray.imzML" -const TEST_IMZML_FILE = "/home/pixel/Documents/Cinvestav_2025/Analisis/Imaging prueba Roya 1/royaimg.imzML" +# const TEST_IMZML_FILE = "/home/pixel/Documents/Cinvestav_2025/Analisis/Imaging prueba Roya 1/royaimg.imzML" # The m/z value to use for creating an image slice. # const MZ_VALUE_FOR_SLICE = 309.06 # BF # const MZ_VALUE_FOR_SLICE = 896.0 # HR2MSI # const MZ_VALUE_FOR_SLICE = 76.03 # I PS -const MZ_VALUE_FOR_SLICE = 473 # ROYA +const MZ_VALUE_FOR_SLICE = 313 # ROYA # const MZ_TOLERANCE = 0.1 # const MZ_TOLERANCE = 1 -const MZ_TOLERANCE = 0.5 +const MZ_TOLERANCE = 0.1 # Coordinates to plot a specific spectrum from imzML const COORDS_TO_PLOT = (50, 50) # Example coordinates (X, Y) @@ -69,8 +71,8 @@ const COORDS_TO_PLOT = (50, 50) # Example coordinates (X, Y) # --- Output Directory --- const RESULTS_DIR = "test/results" -test1 = false -test2 = false +test1 = true +test2 = true test3 = true # =================================================================== @@ -152,7 +154,7 @@ function run_test() try # Get the msi data from the mzml println("Plotting a sample spectrum from $TEST_MZML_FILE...") - msi_data = OpenMSIData(TEST_MZML_FILE) + msi_data = @time OpenMSIData(TEST_MZML_FILE) mz, intensity = GetSpectrum(msi_data, SPECTRUM_TO_PLOT) fig = Figure(size = (800, 600)) @@ -229,10 +231,9 @@ function run_test() if isfile(TEST_IMZML_FILE) # Add spectrum plotting for imzML to match Test Case 1 try - """ # Get the msi data from the imzml println("Plotting a sample spectrum from $TEST_IMZML_FILE...") - msi_data = OpenMSIData(TEST_IMZML_FILE) + msi_data = @time OpenMSIData(TEST_IMZML_FILE) x_coord, y_coord = COORDS_TO_PLOT # Get the x y coordinate spectrum data @@ -275,7 +276,6 @@ function run_test() output_path = joinpath(RESULTS_DIR, "test_imzml_average_spectrum.png") save(output_path, fig) println("SUCCESS: Total spectrum plot saved to $output_path") - """ println("No spectrums tested on this try.") catch e println("ERROR during spectrum plotting in Test Case 3: $e") @@ -284,7 +284,7 @@ function run_test() # Test the plot_slice function try println("Testing plot_slice function on $TEST_IMZML_FILE...") - msi_data = OpenMSIData(TEST_IMZML_FILE) + msi_data = @time OpenMSIData(TEST_IMZML_FILE) @time plot_slice(msi_data, MZ_VALUE_FOR_SLICE, MZ_TOLERANCE, RESULTS_DIR, stage_name="test_imzml_single_slice") # The success message is now inside plot_slice catch e