JuliaMSI
@@ -12,7 +11,7 @@
-
+
@@ -58,7 +57,7 @@
Search for the imzML or mzML file in your system
-
-
+
-
{{msg}}
-
-
@@ -169,25 +168,25 @@
-
-
+
Image visualizer
-
-
+
-
-
{{msgimg}}
@@ -198,17 +197,17 @@TrIQ visualizer
-
-
+
-
-
-
@@ -263,23 +262,23 @@
-
-
{{msgtriq}}
@@ -216,16 +215,16 @@
-
-
-
+
-
+
-
{{msgimg}}
@@ -289,56 +288,56 @@
-
-
+
-
-
-
-
+
-
{{msgtriq}}
-
-
-
+
", line)
+ break
+ end
+ end
+ spectrum_xml = String(take!(spectrum_buffer))
+ break # Found the first spectrum, so we can stop
+ end
+ end
+ finally
+ seek(stream, start_pos) # Always reset stream position
+ end
-"""
- load_imzml_lazy(file_path::String; cache_size=100)
+ if isempty(spectrum_xml)
+ # Fallback based on compression flags if no spectrum tag found
+ return (mz_is_compressed || int_is_compressed) ? :compressed : :uncompressed
+ end
+
+ # Inspect the XML content
+ has_neofx_markers = occursin("encodedLength=\"0\"", spectrum_xml) &&
+ occursin("external encoded length", spectrum_xml)
+
+ has_external_data_markers = occursin("IMS:1000101", spectrum_xml) &&
+ occursin("IMS:1000102", spectrum_xml) &&
+ occursin("IMS:1000103", spectrum_xml)
+
+ if has_neofx_markers
+ return :neofx
+ end
+
+ if mz_is_compressed || int_is_compressed || has_external_data_markers
+ return :compressed
+ end
+
+ return :uncompressed
+end
-Main function to parse an `.imzML`/.ibd file pair and prepare for lazy loading.
-It now returns a unified MSIData object.
-"""
function load_imzml_lazy(file_path::String; cache_size=100)
println("DEBUG: Checking for .imzML file at $file_path")
if !isfile(file_path)
@@ -162,71 +215,64 @@ function load_imzml_lazy(file_path::String; cache_size=100)
try
println("DEBUG: Configuring axes...")
- axis = axes_config_img(stream)
+ param_groups = axes_config_img(stream)
println("DEBUG: Getting image dimensions...")
imgDim = get_img_dimensions(stream)
width, height, num_spectra = imgDim
println("DEBUG: Image dimensions: $(width)x$(height), $num_spectra spectra.")
- mz_config_idx = findfirst(a -> a.Axis == 1, axis)
- 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")
+ # Extract default formats from the parsed param_groups
+ mz_group = nothing
+ int_group = nothing
- # --- NEW PARSING LOGIC based on the old, working code ---
- println("DEBUG: Learning file structure from first spectrum...")
- start_of_spectra_xml = position(stream)
- attr = get_spectrum_attributes(stream, hIbd)
- current_ibd_offset = position(hIbd)
- seek(stream, start_of_spectra_xml)
- println("DEBUG: Initial IBD offset: $current_ibd_offset")
-
- spectra_metadata = Vector{SpectrumMetadata}(undef, num_spectra)
- mz_is_first = attr[3] == 3
-
- println("DEBUG: Parsing metadata for $num_spectra spectra using skip-based method...")
- for k in 1:num_spectra
- # Use skip values learned from the first spectrum, assuming all are identical.
- skip(stream, attr[5]) # Skip to X coordinate value
- val_tag_x = find_tag(stream, r"value=\"(\d+)\"")
- x = parse(Int32, val_tag_x.captures[1])
-
- skip(stream, attr[6]) # Skip to Y coordinate value
- val_tag_y = find_tag(stream, r"value=\"(\d+)\"")
- y = parse(Int32, val_tag_y.captures[1])
-
- skip(stream, attr[7]) # Skip to array length value
- val_tag_len = find_tag(stream, r"value=\"(\d+)\"")
- nPoints = parse(Int32, val_tag_len.captures[1])
-
- mz_len_bytes = nPoints * sizeof(mz_format)
- int_len_bytes = nPoints * sizeof(intensity_format)
-
- local mz_offset, int_offset
- if mz_is_first
- mz_offset = current_ibd_offset
- int_offset = mz_offset + mz_len_bytes
- else
- int_offset = current_ibd_offset
- mz_offset = int_offset + int_len_bytes
+ for group in values(param_groups)
+ if group.Axis == 1
+ mz_group = group
+ elseif group.Axis == 2
+ int_group = group
end
-
- # Create modern SpectrumAsset objects
- 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)
-
- # Advance the offset for the next spectrum's data.
- current_ibd_offset += mz_len_bytes + int_len_bytes
-
- skip(stream, attr[8]) # Skip to the end of the spectrum tag
end
- # --- END OF NEW PARSING LOGIC ---
+
+ if mz_group === nothing || int_group === nothing
+ @warn "Could not find global definitions for m/z and intensity arrays. Using hardcoded defaults (Float64)."
+ default_mz_format = Float64
+ default_intensity_format = Float64
+ mz_is_compressed = false
+ int_is_compressed = false
+ global_mode = UNKNOWN
+ else
+ default_mz_format = mz_group.Format
+ default_intensity_format = int_group.Format
+ mz_is_compressed = mz_group.Packed
+ int_is_compressed = int_group.Packed
+ global_mode = mz_group.Mode != UNKNOWN ? mz_group.Mode : int_group.Mode
+ end
+
+ println("DEBUG: m/z format: $default_mz_format, Intensity format: $default_intensity_format")
+ println("DEBUG: m/z compressed: $mz_is_compressed, Intensity compressed: $int_is_compressed")
+ println("DEBUG: Global mode: $global_mode")
+
+ # --- Parser Selection ---
+ parser_type = determine_parser(stream, mz_is_compressed, int_is_compressed)
+ println("DEBUG: Selected parser: $parser_type")
+
+ local spectra_metadata
+ if parser_type == :neofx
+ println("DEBUG: Using neofx parser.")
+ spectra_metadata = parse_neofx(stream, hIbd, param_groups, width, height, num_spectra,
+ default_mz_format, default_intensity_format,
+ mz_is_compressed, int_is_compressed, global_mode)
+ elseif parser_type == :compressed
+ println("DEBUG: Using compressed parser.")
+ spectra_metadata = parse_compressed(stream, hIbd, param_groups, width, height, num_spectra,
+ default_mz_format, default_intensity_format,
+ mz_is_compressed, int_is_compressed, global_mode)
+ else # :uncompressed
+ println("DEBUG: Using uncompressed parser.")
+ spectra_metadata = parse_uncompressed(stream, hIbd, param_groups, width, height, num_spectra,
+ default_mz_format, default_intensity_format,
+ mz_is_compressed, int_is_compressed, global_mode)
+ end
println("DEBUG: Metadata parsing complete.")
@@ -234,6 +280,9 @@ function load_imzml_lazy(file_path::String; cache_size=100)
println("DEBUG: Building coordinate map...")
coordinate_map = zeros(Int, width, height)
for (idx, meta) in enumerate(spectra_metadata)
+ if idx == 1
+ println("DIAGNOSTIC_WRITE: For index 1, attempting to write to coordinate_map[$(meta.x), $(meta.y)]")
+ end
if 1 <= meta.x <= width && 1 <= meta.y <= height
coordinate_map[meta.x, meta.y] = idx
end
@@ -242,7 +291,7 @@ function load_imzml_lazy(file_path::String; cache_size=100)
close(stream)
- source = ImzMLSource(hIbd, mz_format, intensity_format)
+ 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)
@@ -253,11 +302,342 @@ function load_imzml_lazy(file_path::String; cache_size=100)
end
end
+function parse_uncompressed(stream, hIbd, param_groups, width, height, num_spectra,
+ mz_format, intensity_format, mz_is_compressed, int_is_compressed, global_mode)
+ # Your existing working skip-based parser
+ println("DEBUG: Learning file structure from first spectrum...")
+ start_of_spectra_xml = position(stream)
+ attr = get_spectrum_attributes(stream, hIbd)
+ current_ibd_offset = position(hIbd)
+ seek(stream, start_of_spectra_xml)
+ println("DEBUG: Initial IBD offset: $current_ibd_offset")
+
+ spectra_metadata = Vector{SpectrumMetadata}(undef, num_spectra)
+ mz_is_first = attr[3] == 3
+
+ for k in 1:num_spectra
+ # Store the start position of this spectrum for mode detection
+ spectrum_start_pos = position(stream)
+
+ # Use skip values learned from the first spectrum
+ skip(stream, attr[5]) # Skip to X coordinate value
+ val_tag_x = find_tag(stream, r"value=\"(\d+)\"")
+ x = parse(Int32, val_tag_x.captures[1])
+
+ skip(stream, attr[6]) # Skip to Y coordinate value
+ val_tag_y = find_tag(stream, r"value=\"(\d+)\"")
+ y = parse(Int32, val_tag_y.captures[1])
+
+ skip(stream, attr[7]) # Skip to array length value
+ val_tag_len = find_tag(stream, r"value=\"(\d+)\"")
+ nPoints = parse(Int32, val_tag_len.captures[1])
+
+ # For uncompressed data, use simple calculation
+ mz_len_bytes = nPoints * sizeof(mz_format)
+ int_len_bytes = nPoints * sizeof(intensity_format)
+
+ local mz_offset, int_offset
+ if mz_is_first
+ mz_offset = current_ibd_offset
+ int_offset = mz_offset + mz_len_bytes
+ else
+ int_offset = current_ibd_offset
+ mz_offset = int_offset + int_len_bytes
+ end
+
+ # Mode detection from spectrum XML
+ current_pos = position(stream)
+ seek(stream, spectrum_start_pos)
+ spectrum_buffer = IOBuffer()
+ line = ""
+ while !eof(stream)
+ line = readline(stream)
+ write(spectrum_buffer, line)
+ if occursin("", line)
+ break
+ end
+ end
+ spectrum_xml = String(take!(spectrum_buffer))
+
+ spectrum_mode = global_mode
+ if occursin("MS:1000127", spectrum_xml)
+ spectrum_mode = CENTROID
+ elseif occursin("MS:1000128", spectrum_xml)
+ spectrum_mode = PROFILE
+ end
+ seek(stream, current_pos)
+
+ # Create SpectrumAsset objects
+ 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, "", spectrum_mode, mz_asset, int_asset)
+
+ current_ibd_offset += mz_len_bytes + int_len_bytes
+ skip(stream, attr[8]) # Skip to the end of the spectrum tag
+ end
+
+ return spectra_metadata
+end
+
+function parse_compressed(stream, hIbd, param_groups, width, height, num_spectra,
+ default_mz_format, default_intensity_format,
+ mz_is_compressed, int_is_compressed, global_mode)
+ # New parser for compressed data
+ spectra_metadata = Vector{SpectrumMetadata}(undef, num_spectra)
+
+ for k in 1:num_spectra
+ # Read the full spectrum XML block
+ spectrum_buffer = IOBuffer()
+ line = ""
+ while !eof(stream)
+ line = readline(stream)
+ if occursin("", 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 = []
+
+ # Find all binaryDataArray blocks
+ array_matches = eachmatch(r""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 = 0
+ if array_len_cv_match !== nothing
+ array_length = parse(Int32, array_len_cv_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])
+ end
+ end
+
+ # Get encoded_length
+ encoded_len_cv_match = match(r"IMS:1000104.*?value=\"(\d+)\"", array_xml)
+ encoded_length = 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 = 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
+ 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)
+
+ if length(mz_data) != 1 || length(int_data) != 1
+ error("Spectrum $k: Expected exactly one m/z and one intensity array")
+ end
+
+ 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,
+ int_is_compressed ? int_info.encoded_length : int_info.array_length, :intensity)
+
+ spectra_metadata[k] = SpectrumMetadata(x, y, "", spectrum_mode, mz_asset, int_asset)
+ end
+
+ return spectra_metadata
+end
+
+function parse_neofx(stream, hIbd, param_groups, width, height, num_spectra,
+ default_mz_format, default_intensity_format,
+ mz_is_compressed, int_is_compressed, global_mode)
+ # New parser for compressed data
+ spectra_metadata = Vector{SpectrumMetadata}(undef, num_spectra)
+
+ for k in 1:num_spectra
+ # Read the full spectrum XML block
+ spectrum_buffer = IOBuffer()
+ line = ""
+ while !eof(stream)
+ line = readline(stream)
+ if occursin("", 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 = []
+
+ # Find all binaryDataArray blocks
+ array_matches = eachmatch(r""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 = 0
+ if array_len_cv_match !== nothing
+ array_length = parse(Int32, array_len_cv_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])
+ end
+ end
+
+ # Get encoded_length
+ encoded_len_cv_match = match(r"IMS:1000104.*?value=\"(\d+)\"", array_xml)
+ encoded_length = 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 = 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
+ 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)
+
+ if length(mz_data) != 1 || length(int_data) != 1
+ error("Spectrum $k: Expected exactly one m/z and one intensity array")
+ end
+
+ 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,
+ int_is_compressed ? int_info.encoded_length : int_info.array_length, :intensity)
+
+ spectra_metadata[k] = SpectrumMetadata(x, y, "", spectrum_mode, mz_asset, int_asset)
+ end
+
+ return spectra_metadata
+end
+
# --- End of content from imzML.jl ---
-
-# --- Start of content from Imaging_Normalization.jl ---
-
# =============================================================================
#
# Image Slice Extraction
@@ -446,7 +826,8 @@ function plot_slice(msi_data::MSIData, mass::Real, tolerance::Real, output_dir::
fig = Figure(size = (600, 500))
ax = CairoMakie.Axis(fig[1, 1],
aspect=DataAspect(),
- title=@sprintf("Slice for m/z: %.2f", mass)
+ title=@sprintf("Slice for m/z: %.2f", mass),
+ yreversed=true
)
hidedecorations!(ax)
@@ -985,4 +1366,4 @@ function generate_palette(colorscheme, n_colors=256)
return palette
end
-const ViridisPalette = generate_palette(ColorSchemes.viridis)
+const ViridisPalette = generate_palette(ColorSchemes.viridis)
\ No newline at end of file
diff --git a/test/run_tests.jl b/test/run_tests.jl
index d8c9de9..9c7d829 100644
--- a/test/run_tests.jl
+++ b/test/run_tests.jl
@@ -52,15 +52,19 @@ 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_paper_spray/Imaging_paper_spray.imzML" #profile
# const TEST_IMZML_FILE = "/home/pixel/Documents/Cinvestav_2025/Analisis/Imaging prueba Roya 1/royaimg.imzML"
+# const TEST_IMZML_FILE = "/home/pixel/Documents/Cinvestav_2025/Analisis/salida/ltpmsi-chilli.imzML" # centroid aparently?
+# const TEST_IMZML_FILE = "/home/pixel/Documents/Cinvestav_2025/Analisis/salida/Stomach_DHB_compressed.imzML" # centroid compressed
+const TEST_IMZML_FILE = "/home/pixel/Documents/Cinvestav_2025/Analisis/salida/Stomach_DHB_uncompressed.imzML" # centroid
# 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 = 313 # ROYA
+# const MZ_VALUE_FOR_SLICE = 313 # ROYA
+const MZ_VALUE_FOR_SLICE = 100 # advanced processing
# const MZ_TOLERANCE = 0.1
# const MZ_TOLERANCE = 1
const MZ_TOLERANCE = 0.1
@@ -71,8 +75,8 @@ const COORDS_TO_PLOT = (50, 50) # Example coordinates (X, Y)
# --- Output Directory ---
const RESULTS_DIR = "test/results"
-test1 = true
-test2 = true
+test1 = false
+test2 = false
test3 = true
# ===================================================================
@@ -135,6 +139,30 @@ function validate_msi_data(filepath::String)
end
end
+function debug_xml_parsing(file_path::String)
+ println("=== DEBUG XML PARSING ===")
+ stream = open(file_path, "r")
+
+ # Find and print the first spectrum
+ while !eof(stream)
+ line = readline(stream)
+ if occursin("", line)
+ line = readline(stream)
+ spectrum_xml *= line
+ end
+ println("SPECTRUM XML:")
+ println(spectrum_xml)
+ break
+ end
+ end
+ close(stream)
+ println("=== END DEBUG ===")
+end
+
# ===================================================================
# TEST RUNNER
@@ -229,16 +257,24 @@ function run_test()
# Also run tests for plotting spectrum and image slice
if isfile(TEST_IMZML_FILE)
+ debug_xml_parsing(TEST_IMZML_FILE)
+ msi_data = @time OpenMSIData(TEST_IMZML_FILE)
+ precompute_analytics(msi_data)
# 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 = @time OpenMSIData(TEST_IMZML_FILE)
- x_coord, y_coord = COORDS_TO_PLOT
+
+ # Use coordinates from the first spectrum in the metadata
+ first_spectrum_meta = msi_data.spectra_metadata[1]
+ x_coord = first_spectrum_meta.x
+ y_coord = first_spectrum_meta.y
+
+ println("DIAGNOSTIC_READ: Reading from coordinate_map[$x_coord, $y_coord]. Value is $(msi_data.coordinate_map[x_coord, y_coord]).")
# Get the x y coordinate spectrum data
- mz, intensity = GetSpectrum(msi_data, x_coord, y_coord)
-
+ mz, intensity = GetSpectrum(msi_data, Int(x_coord), Int(y_coord))
# Plot the data
fig = Figure(size = (800, 600))
ax = Axis(fig[1, 1], xlabel="m/z", ylabel="Intensity", title="Spectrum at ($x_coord, $y_coord) from $(basename(TEST_IMZML_FILE))")
@@ -276,7 +312,7 @@ 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.")
+ # println("No spectrums tested on this try.")
catch e
println("ERROR during spectrum plotting in Test Case 3: $e")
end
@@ -284,8 +320,13 @@ function run_test()
# Test the plot_slice function
try
println("Testing plot_slice function on $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")
+
+ # Use the base peak m/z from the first spectrum for the slice
+ # Or a more general approach: use the global max intensity m/z
+ # For now, let's use the base peak m/z of the first spectrum
+ slice_mz_value = msi_data.spectrum_stats_df.BasePeakMZ[1]
+
+ @time plot_slice(msi_data, slice_mz_value, MZ_TOLERANCE, RESULTS_DIR, stage_name="test_imzml_single_slice")
# The success message is now inside plot_slice
catch e
println("ERROR during plot_slice test in Test Case 3: $e")
-
-
{{msgimgComp}}
@@ -348,34 +347,34 @@
-
-
+
-
-
-
@@ -411,6 +410,4 @@
-
-
{{msgtriqComp}}