diff --git a/Manifest.toml b/Manifest.toml index 4dfa48d..4ec7d85 100644 --- a/Manifest.toml +++ b/Manifest.toml @@ -2,7 +2,7 @@ julia_version = "1.11.7" manifest_format = "2.0" -project_hash = "9b36a3561cfc9927071202de68baf8c2b0f02a5c" +project_hash = "f5ab935406d0b7b7e657445ecfa45944ea71f247" [[deps.ATK_jll]] deps = ["Artifacts", "Glib_jll", "JLLWrappers", "Libdl"] diff --git a/Project.toml b/Project.toml index 816068c..5d7b777 100644 --- a/Project.toml +++ b/Project.toml @@ -32,8 +32,10 @@ NativeFileDialog = "e1fe445b-aa65-4df4-81c1-2041507f0fd4" NaturalSort = "c020b1a1-e9b0-503a-9c33-f039bfc54a85" PlotlyBase = "a03496cd-edff-5a9b-9e67-9cda94a718b5" Printf = "de0858da-6303-5e67-8744-51eddeeeb8d7" +ProgressMeter = "92933f4c-e287-5a05-a399-4b506db050ca" SavitzkyGolay = "c4bf5708-b6a6-4fbe-bcd0-6850ed671584" Serialization = "9e88b42a-f829-5b0c-bbe9-9e923198166b" Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2" StatsBase = "2913bbd2-ae8a-5f71-8c99-4fb6c76f3a91" StipplePlotly = "ec984513-233d-481d-95b0-a3b58b97af2b" +UUIDs = "cf7118a7-6976-5b1a-9a39-7adc72f591a4" diff --git a/README.md b/README.md index 2affa30..489bd0f 100644 --- a/README.md +++ b/README.md @@ -22,7 +22,7 @@ https://codeberg.org/LabABI/JuliaMSI ~/Downloads/JuliaMSI-main/juliamsi 2. Without entering the Julia environment, launch the project in your terminal with the following command (which works for all operating systems): ``` - julia --project=. start_MSI_GUI.jl + julia --threads auto --project=. start_MSI_GUI.jl ``` 3. After the script has finished loading, you can open a [page](http://127.0.0.1:1481/) in your browser with the web app running. diff --git a/app.jl b/app.jl index 1d6ff84..1e96ee7 100644 --- a/app.jl +++ b/app.jl @@ -213,8 +213,8 @@ end # == Batch Processing & Registry Variables == @private registry_init_done = false @in selected_files = String[] - @out available_folders = String[] - @out image_available_folders = String[] + @in available_folders = String[] + @in image_available_folders = String[] @out registry_path = abspath(joinpath(@__DIR__, "public", "registry.json")) # Progress reporting @out overall_progress = 0.0 @@ -317,7 +317,7 @@ end margin=attr(l=0,r=0,t=120,b=0,pad=0) ) # Dummy 2D scatter plot - traceSpectra=PlotlyBase.stem(x=Vector{Float64}(), y=Vector{Float64}(),marker=attr(size=1, color="blue", opacity=0.1)) + traceSpectra=PlotlyBase.scatter(x=Vector{Float64}(), y=Vector{Float64}(), mode="lines", marker=attr(size=1, color="blue", opacity=0.1)) # Create conection to frontend @out plotdata=[traceSpectra] @out plotlayout=layoutSpectra @@ -410,113 +410,111 @@ end progress = true msg = "Opening file: $(basename(picked_route))..." - @async begin - try - dataset_name = replace(basename(picked_route), r"(\.(imzML|imzml|mzML|mzml))$"i => "") - registry = load_registry(registry_path) - existing_entry = get(registry, dataset_name, nothing) + try + dataset_name = replace(basename(picked_route), r"(\.(imzML|imzml|mzML|mzml))$"i => "") + registry = load_registry(registry_path) + existing_entry = get(registry, dataset_name, nothing) - # --- Fast Load Path --- - is_same_file = (existing_entry !== nothing && existing_entry["source_path"] == picked_route) - if is_same_file && !isempty(get(existing_entry, "metadata", Dict())) - msg = "Fast loading pre-processed file: $(dataset_name)" - println(msg) - - full_route = existing_entry["source_path"] - metadata_rows = existing_entry["metadata"]["summary"] - - dims_str = first(filter(r -> r["parameter"] == "Image Dimensions", metadata_rows))["value"] - dims = parse.(Int, split(dims_str, " x ")) - imgWidth, imgHeight = dims[1], dims[2] - - msi_data = nothing # Ensure data is not held in memory - log_memory_usage("Fast Load (msi_data cleared)", msi_data) - btnMetadataDisable = false - btnStartDisable = false - btnPlotDisable = false - btnSpectraDisable = false - SpectraEnabled = true - selected_folder_main = dataset_name - - # Update folder lists in UI - all_folders = sort(collect(keys(registry)), lt=natural) - img_folders = filter(folder -> get(get(registry, folder, Dict()), "is_imzML", false), all_folders) - available_folders = deepcopy(all_folders) - image_available_folders = deepcopy(img_folders) - - msg = "Successfully loaded pre-processed dataset: $(dataset_name)" - progress = false - return - end - - # --- Full Load Path --- - msg = "Performing first-time analysis for: $(basename(picked_route))..." - local local_full_route - if endswith(picked_route, r"imzml"i) - local_full_route = replace(picked_route, r"\.imzml$"i => ".imzML") - if picked_route != local_full_route - mv(picked_route, local_full_route, force=true) - end - else - local_full_route = picked_route - end - full_route = local_full_route - - sTime = time() - loaded_data = OpenMSIData(local_full_route) - is_imzML = loaded_data.source isa ImzMLSource + # --- Fast Load Path --- + is_same_file = (existing_entry !== nothing && existing_entry["source_path"] == picked_route) + if is_same_file && !isempty(get(existing_entry, "metadata", Dict())) + msg = "Fast loading pre-processed file: $(dataset_name)" + println(msg) - precompute_analytics(loaded_data) + full_route = existing_entry["source_path"] + metadata_rows = existing_entry["metadata"]["summary"] - metadata_columns = [ - Dict("name" => "parameter", "label" => "Parameter", "field" => "parameter", "align" => "left"), - Dict("name" => "value", "label" => "Value", "field" => "value", "align" => "left"), - ] - summary_stats = extract_metadata(loaded_data, local_full_route) - metadata_rows = summary_stats["summary"] - btnMetadataDisable = isempty(metadata_rows) + dims_str = first(filter(r -> r["parameter"] == "Image Dimensions", metadata_rows))["value"] + dims = parse.(Int, split(dims_str, " x ")) + imgWidth, imgHeight = dims[1], dims[2] - w, h = loaded_data.image_dims - imgWidth, imgHeight = w > 0 ? (w, h) : (500, 500) - - update_registry(registry_path, dataset_name, local_full_route, summary_stats, is_imzML) + msi_data = nothing # Ensure data is not held in memory + log_memory_usage("Fast Load (msi_data cleared)", msi_data) + btnMetadataDisable = false + btnStartDisable = false + btnPlotDisable = false + btnSpectraDisable = false + SpectraEnabled = true + selected_folder_main = dataset_name # Update folder lists in UI - registry = load_registry(registry_path) all_folders = sort(collect(keys(registry)), lt=natural) img_folders = filter(folder -> get(get(registry, folder, Dict()), "is_imzML", false), all_folders) available_folders = deepcopy(all_folders) image_available_folders = deepcopy(img_folders) - selected_folder_main = dataset_name - msi_data = loaded_data - log_memory_usage("Full Load", msi_data) - - eTime = round(time() - sTime, digits=3) - msg = "Active file loaded in $(eTime) seconds. Dataset '$(dataset_name)' is ready for analysis." - - btnStartDisable = false - btnPlotDisable = false - btnSpectraDisable = false - SpectraEnabled = true - - catch e - msi_data = nothing - msg = "Error loading active file: $e" - warning_msg = true - btnStartDisable = true - btnSpectraDisable = true - SpectraEnabled = false - btnMetadataDisable = true - @error "File loading failed" exception=(e, catch_backtrace()) - finally - GC.gc() - if Sys.islinux() - ccall(:malloc_trim, Int32, (Int32,), 0) - end + msg = "Successfully loaded pre-processed dataset: $(dataset_name)" progress = false - progressSpectraPlot = false + return end + + # --- Full Load Path --- + msg = "Performing first-time analysis for: $(basename(picked_route))..." + local local_full_route + if endswith(picked_route, r"imzml"i) + local_full_route = replace(picked_route, r"\.imzml$"i => ".imzML") + if picked_route != local_full_route + mv(picked_route, local_full_route, force=true) + end + else + local_full_route = picked_route + end + full_route = local_full_route + + sTime = time() + loaded_data = OpenMSIData(local_full_route) + is_imzML = loaded_data.source isa ImzMLSource + + precompute_analytics(loaded_data) + + metadata_columns = [ + Dict("name" => "parameter", "label" => "Parameter", "field" => "parameter", "align" => "left"), + Dict("name" => "value", "label" => "Value", "field" => "value", "align" => "left"), + ] + summary_stats = extract_metadata(loaded_data, local_full_route) + metadata_rows = summary_stats["summary"] + btnMetadataDisable = isempty(metadata_rows) + + w, h = loaded_data.image_dims + imgWidth, imgHeight = w > 0 ? (w, h) : (500, 500) + + update_registry(registry_path, dataset_name, local_full_route, summary_stats, is_imzML) + + # Update folder lists in UI + registry = load_registry(registry_path) + all_folders = sort(collect(keys(registry)), lt=natural) + img_folders = filter(folder -> get(get(registry, folder, Dict()), "is_imzML", false), all_folders) + available_folders = deepcopy(all_folders) + image_available_folders = deepcopy(img_folders) + + selected_folder_main = dataset_name + msi_data = loaded_data + log_memory_usage("Full Load", msi_data) + + eTime = round(time() - sTime, digits=3) + msg = "Active file loaded in $(eTime) seconds. Dataset '$(dataset_name)' is ready for analysis." + + btnStartDisable = false + btnPlotDisable = false + btnSpectraDisable = false + SpectraEnabled = true + + catch e + msi_data = nothing + msg = "Error loading active file: $e" + warning_msg = true + btnStartDisable = true + btnSpectraDisable = true + SpectraEnabled = false + btnMetadataDisable = true + @error "File loading failed" exception=(e, catch_backtrace()) + finally + GC.gc() + if Sys.islinux() + ccall(:malloc_trim, Int32, (Int32,), 0) + end + progress = false + progressSpectraPlot = false end end @@ -620,34 +618,32 @@ end btnConvertDisable = true msg_conversion = "Starting conversion process..." - @async begin - try - sTime = time() - target_imzml = replace(mzml_full_route, r"\.(mzml|mzML)$" => ".imzML") - - msg_conversion = "Converting $(basename(mzml_full_route)) to $(basename(target_imzml))... This may take a while." + try + sTime = time() + target_imzml = replace(mzml_full_route, r"\.(mzml|mzML)$" => ".imzML") + + msg_conversion = "Converting $(basename(mzml_full_route)) to $(basename(target_imzml))... This may take a while." - success = ImportMzmlFile(mzml_full_route, sync_full_route, target_imzml) + success = ImportMzmlFile(mzml_full_route, sync_full_route, target_imzml) - fTime = time() - eTime = round(fTime - sTime, digits=3) + fTime = time() + eTime = round(fTime - sTime, digits=3) - if success - msg_conversion = "Conversion successful in $(eTime) seconds. Output file: $(basename(target_imzml))" - else - msg_conversion = "Conversion failed after $(eTime) seconds. Check console for errors." - warning_msg = true - end - - catch e - msg_conversion = "An error occurred during conversion: $e" + if success + msg_conversion = "Conversion successful in $(eTime) seconds. Output file: $(basename(target_imzml))" + else + msg_conversion = "Conversion failed after $(eTime) seconds. Check console for errors." warning_msg = true - @error "Conversion failed" exception=(e, catch_backtrace()) - finally - progress_conversion = false - # Re-enable button if files are still selected - btnConvertDisable = isempty(mzml_full_route) || isempty(sync_full_route) end + + catch e + msg_conversion = "An error occurred during conversion: $e" + warning_msg = true + @error "Conversion failed" exception=(e, catch_backtrace()) + finally + progress_conversion = false + # Re-enable button if files are still selected + btnConvertDisable = isempty(mzml_full_route) || isempty(sync_full_route) end end @@ -672,174 +668,171 @@ end current_registry_path = registry_path println("starting main process with $(length(current_selected_files)) files") - - @async begin - total_time_start = time() - try - # --- 1. Parameter Validation --- - if isempty(current_selected_files) - progress_message = "No .imzML files in batch. Please add files first." - warning_msg = true - println(progress_message) - return - end - - masses = Float64[] - try - masses = [parse(Float64, strip(m)) for m in split(current_nmass, ',', keepempty=false)] - catch e - progress_message = "Invalid m/z value(s). Please provide a comma-separated list of numbers. Error: $e" - warning_msg = true - return - end - - if isempty(masses) - progress_message = "No valid m/z values found. Please provide comma-separated positive numbers." - warning_msg = true - return - end - - # --- 2. Batch Processing Loop --- - num_files = length(current_selected_files) - total_steps = num_files - current_step = 0 - errors = Dict("load_errors" => String[], "slice_errors" => String[], "io_errors" => String[]) - newly_created_folders = String[] - files_without_mask = 0 - - for (file_idx, file_path) in enumerate(current_selected_files) - progress_message = "Processing file $file_idx/$num_files: $(basename(file_path))" - overall_progress = current_step / total_steps - - all_params = ( - tolerance = current_tol, - colorL = current_color_level, - triqE = current_triq_enabled, - triqP = current_triq_prob, - medianF = current_mfilter_enabled, - registry = current_registry_path, - fileIdx = file_idx, - nFiles = num_files - ) - - success, error_msg = process_file_safely(file_path, masses, all_params, progress_message, overall_progress, use_mask=current_mask_enabled) - - if !success - push!(errors["load_errors"], error_msg) - else - push!(newly_created_folders, replace(basename(file_path), r"\.imzML$"i => "")) - end - current_step += 1 - end - - # --- 3. Final Report --- - total_time_end = round(time() - total_time_start, digits=3) - - registry = load_registry(current_registry_path) - all_folders = sort(collect(keys(registry)), lt=natural) - img_folders = filter(folder -> get(get(registry, folder, Dict()), "is_imzML", false), all_folders) - available_folders = deepcopy(all_folders) - image_available_folders = deepcopy(img_folders) - - if !isempty(newly_created_folders) - selected_folder_main = first(newly_created_folders) - end - - successful_files = length(newly_created_folders) - total_errors = sum(length, values(errors)) - - if total_errors == 0 - msg = "Successfully processed all $(successful_files) file(s) in $(total_time_end) seconds." - else - msg = "Batch completed in $(total_time_end) seconds with $(total_errors) error(s)." - warning_msg = true - end - - mask_summary = current_mask_enabled ? "\nFiles processed without a mask: $(files_without_mask)" : "" - - batch_summary = """ - Processed $(successful_files)/$(num_files) files successfully. - $(mask_summary) - - Errors by category: - • Load failures: $(length(errors["load_errors"])) - • Slice generation: $(length(errors["slice_errors"])) - • I/O issues: $(length(errors["io_errors"])) - - Detailed errors: - $(join(vcat(values(errors)...), "\n")) - """ - showBatchSummary = true - - # Update UI to display the last generated image - if !isempty(newly_created_folders) - timestamp = string(time_ns()) - folder_path = joinpath("public", selected_folder_main) - - if current_triq_enabled - triq_files = filter(filename -> startswith(filename, "TrIQ_") && endswith(filename, ".bmp"), readdir(folder_path)) - col_triq_files = filter(filename -> startswith(filename, "colorbar_TrIQ_") && endswith(filename, ".png"), readdir(folder_path)) - - if !isempty(triq_files) - latest_triq = triq_files[argmax([mtime(joinpath(folder_path, f)) for f in triq_files])] - current_triq = latest_triq - imgIntT = "/$(selected_folder_main)/$(current_triq)?t=$(timestamp)" - plotdataImgT, plotlayoutImgT, _, _ = loadImgPlot(imgIntT) - text_nmass = replace(current_triq, r"TrIQ_|.bmp" => "") - msgtriq = "TrIQ m/z: $(replace(text_nmass, "_" => "."))" - - if !isempty(col_triq_files) - latest_col_triq = col_triq_files[argmax([mtime(joinpath(folder_path, f)) for f in col_triq_files])] - current_col_triq = latest_col_triq - colorbarT = "/$(selected_folder_main)/$(current_col_triq)?t=$(timestamp)" - else - colorbarT = "" - end - selectedTab = "tab1" - end - else # Not TrIQ enabled, display regular MSI image - msi_files = filter(filename -> startswith(filename, "MSI_") && endswith(filename, ".bmp"), readdir(folder_path)) - col_msi_files = filter(filename -> startswith(filename, "colorbar_MSI_") && endswith(filename, ".png"), readdir(folder_path)) - - if !isempty(msi_files) - latest_msi = msi_files[argmax([mtime(joinpath(folder_path, f)) for f in msi_files])] - current_msi = latest_msi - imgInt = "/$(selected_folder_main)/$(current_msi)?t=$(timestamp)" - plotdataImg, plotlayoutImg, _, _ = loadImgPlot(imgInt) - text_nmass = replace(current_msi, r"MSI_|.bmp" => "") - msgimg = "m/z: $(replace(text_nmass, "_" => "."))" - - if !isempty(col_msi_files) - latest_col_msi = col_msi_files[argmax([mtime(joinpath(folder_path, f)) for f in col_msi_files])] - current_col_msi = latest_col_msi - colorbar = "/$(selected_folder_main)/$(current_col_msi)?t=$(timestamp)" - else - colorbar = "" - end - selectedTab = "tab0" - end - end - end - - catch e - println("Error in main process: $e") - msg = "Batch processing failed: $e" + total_time_start = time() + try + # --- 1. Parameter Validation --- + if isempty(current_selected_files) + progress_message = "No .imzML files in batch. Please add files first." warning_msg = true - @error "Main process failed" exception=(e, catch_backtrace()) - finally - # --- UI State Reset --- - progress = false - btnStartDisable = false - btnPlotDisable = false - btnOpticalDisable = false - btnSpectraDisable = false - SpectraEnabled = true - overall_progress = 0.0 - println("Done") - GC.gc() - if Sys.islinux() - ccall(:malloc_trim, Int32, (Int32,), 0) + println(progress_message) + return + end + + masses = Float64[] + try + masses = [parse(Float64, strip(m)) for m in split(current_nmass, ',', keepempty=false)] + catch e + progress_message = "Invalid m/z value(s). Please provide a comma-separated list of numbers. Error: $e" + warning_msg = true + return + end + + if isempty(masses) + progress_message = "No valid m/z values found. Please provide comma-separated positive numbers." + warning_msg = true + return + end + + # --- 2. Batch Processing Loop --- + num_files = length(current_selected_files) + total_steps = num_files + current_step = 0 + errors = Dict("load_errors" => String[], "slice_errors" => String[], "io_errors" => String[]) + newly_created_folders = String[] + files_without_mask = 0 + + for (file_idx, file_path) in enumerate(current_selected_files) + progress_message = "Processing file $file_idx/$num_files: $(basename(file_path))" + overall_progress = current_step / total_steps + + all_params = ( + tolerance = current_tol, + colorL = current_color_level, + triqE = current_triq_enabled, + triqP = current_triq_prob, + medianF = current_mfilter_enabled, + registry = current_registry_path, + fileIdx = file_idx, + nFiles = num_files + ) + + success, error_msg = process_file_safely(file_path, masses, all_params, progress_message, overall_progress, use_mask=current_mask_enabled) + + if !success + push!(errors["load_errors"], error_msg) + else + push!(newly_created_folders, replace(basename(file_path), r"\.imzML$"i => "")) end + current_step += 1 + end + + # --- 3. Final Report --- + total_time_end = round(time() - total_time_start, digits=3) + + registry = load_registry(current_registry_path) + all_folders = sort(collect(keys(registry)), lt=natural) + img_folders = filter(folder -> get(get(registry, folder, Dict()), "is_imzML", false), all_folders) + available_folders = deepcopy(all_folders) + image_available_folders = deepcopy(img_folders) + + if !isempty(newly_created_folders) + selected_folder_main = first(newly_created_folders) + end + + successful_files = length(newly_created_folders) + total_errors = sum(length, values(errors)) + + if total_errors == 0 + msg = "Successfully processed all $(successful_files) file(s) in $(total_time_end) seconds." + else + msg = "Batch completed in $(total_time_end) seconds with $(total_errors) error(s)." + warning_msg = true + end + + mask_summary = current_mask_enabled ? "\nFiles processed without a mask: $(files_without_mask)" : "" + + batch_summary = """ + Processed $(successful_files)/$(num_files) files successfully. + $(mask_summary) + + Errors by category: + • Load failures: $(length(errors["load_errors"])) + • Slice generation: $(length(errors["slice_errors"])) + • I/O issues: $(length(errors["io_errors"])) + + Detailed errors: + $(join(vcat(values(errors)...), "\n")) + """ + showBatchSummary = true + + # Update UI to display the last generated image + if !isempty(newly_created_folders) + timestamp = string(time_ns()) + folder_path = joinpath("public", selected_folder_main) + + if current_triq_enabled + triq_files = filter(filename -> startswith(filename, "TrIQ_") && endswith(filename, ".bmp"), readdir(folder_path)) + col_triq_files = filter(filename -> startswith(filename, "colorbar_TrIQ_") && endswith(filename, ".png"), readdir(folder_path)) + + if !isempty(triq_files) + latest_triq = triq_files[argmax([mtime(joinpath(folder_path, f)) for f in triq_files])] + current_triq = latest_triq + imgIntT = "/$(selected_folder_main)/$(current_triq)?t=$(timestamp)" + plotdataImgT, plotlayoutImgT, _, _ = loadImgPlot(imgIntT) + text_nmass = replace(current_triq, r"TrIQ_|.bmp" => "") + msgtriq = "TrIQ m/z: $(replace(text_nmass, "_" => "."))" + + if !isempty(col_triq_files) + latest_col_triq = col_triq_files[argmax([mtime(joinpath(folder_path, f)) for f in col_triq_files])] + current_col_triq = latest_col_triq + colorbarT = "/$(selected_folder_main)/$(current_col_triq)?t=$(timestamp)" + else + colorbarT = "" + end + selectedTab = "tab1" + end + else # Not TrIQ enabled, display regular MSI image + msi_files = filter(filename -> startswith(filename, "MSI_") && endswith(filename, ".bmp"), readdir(folder_path)) + col_msi_files = filter(filename -> startswith(filename, "colorbar_MSI_") && endswith(filename, ".png"), readdir(folder_path)) + + if !isempty(msi_files) + latest_msi = msi_files[argmax([mtime(joinpath(folder_path, f)) for f in msi_files])] + current_msi = latest_msi + imgInt = "/$(selected_folder_main)/$(current_msi)?t=$(timestamp)" + plotdataImg, plotlayoutImg, _, _ = loadImgPlot(imgInt) + text_nmass = replace(current_msi, r"MSI_|.bmp" => "") + msgimg = "m/z: $(replace(text_nmass, "_" => "."))" + + if !isempty(col_msi_files) + latest_col_msi = col_msi_files[argmax([mtime(joinpath(folder_path, f)) for f in col_msi_files])] + current_col_msi = latest_col_msi + colorbar = "/$(selected_folder_main)/$(current_col_msi)?t=$(timestamp)" + else + colorbar = "" + end + selectedTab = "tab0" + end + end + end + + catch e + println("Error in main process: $e") + msg = "Batch processing failed: $e" + warning_msg = true + @error "Main process failed" exception=(e, catch_backtrace()) + finally + # --- UI State Reset --- + progress = false + btnStartDisable = false + btnPlotDisable = false + btnOpticalDisable = false + btnSpectraDisable = false + SpectraEnabled = true + overall_progress = 0.0 + println("Done") + GC.gc() + if Sys.islinux() + ccall(:malloc_trim, Int32, (Int32,), 0) end end end @@ -856,60 +849,58 @@ end btnStartDisable = true msg = "Loading plot for $(selected_folder_main)..." - @async begin - 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 - - local mask_path_for_plot::Union{String, Nothing} = nothing - if maskEnabled && get(entry, "has_mask", false) - mask_path_for_plot = get(entry, "mask_path", "") - if !isfile(mask_path_for_plot) - @warn "Mask not found for plotting: $(mask_path_for_plot). Plotting without mask." - mask_path_for_plot = nothing - end - end - - plotdata, plotlayout, xSpectraMz, ySpectraMz = meanSpectrumPlot(msi_data, selected_folder_main, mask_path=mask_path_for_plot) - selectedTab = "tab2" - fTime = time() - eTime = round(fTime - sTime, digits=3) - msg = "Plot loaded in $(eTime) seconds" - log_memory_usage("Mean Plot Generated", msi_data) - catch e - msg = "Could not generate mean spectrum plot: $e" + if target_path == "unknown (manually added)" + msg = "Dataset selected contained no route." warning_msg = true - @error "Mean spectrum plotting failed" exception=(e, catch_backtrace()) - finally - progressSpectraPlot = false - btnPlotDisable = false - btnSpectraDisable = false - btnStartDisable = false - GC.gc() - if Sys.islinux() - ccall(:malloc_trim, Int32, (Int32,), 0) + 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 + + local mask_path_for_plot::Union{String, Nothing} = nothing + if maskEnabled && get(entry, "has_mask", false) + mask_path_for_plot = get(entry, "mask_path", "") + if !isfile(mask_path_for_plot) + @warn "Mask not found for plotting: $(mask_path_for_plot). Plotting without mask." + mask_path_for_plot = nothing + end + end + + plotdata, plotlayout, xSpectraMz, ySpectraMz = meanSpectrumPlot(msi_data, selected_folder_main, mask_path=mask_path_for_plot) + selectedTab = "tab2" + fTime = time() + eTime = round(fTime - sTime, digits=3) + msg = "Plot loaded in $(eTime) seconds" + log_memory_usage("Mean Plot Generated", msi_data) + catch e + msg = "Could not generate mean spectrum plot: $e" + warning_msg = true + @error "Mean spectrum plotting failed" exception=(e, catch_backtrace()) + finally + progressSpectraPlot = false + btnPlotDisable = false + btnSpectraDisable = false + btnStartDisable = false + GC.gc() + if Sys.islinux() + ccall(:malloc_trim, Int32, (Int32,), 0) + end end end @@ -925,60 +916,58 @@ end btnStartDisable = true msg = "Loading total spectrum plot for $(selected_folder_main)..." - @async begin - 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 - - local mask_path_for_plot::Union{String, Nothing} = nothing - if maskEnabled && get(entry, "has_mask", false) - mask_path_for_plot = get(entry, "mask_path", "") - if !isfile(mask_path_for_plot) - @warn "Mask not found for plotting: $(mask_path_for_plot). Plotting without mask." - mask_path_for_plot = nothing - end - end - - plotdata, plotlayout, xSpectraMz, ySpectraMz = sumSpectrumPlot(msi_data, selected_folder_main, mask_path=mask_path_for_plot) - selectedTab = "tab2" - fTime = time() - eTime = round(fTime - sTime, digits=3) - msg = "Total plot loaded in $(eTime) seconds" - log_memory_usage("Sum Plot Generated", msi_data) - catch e - msg = "Could not generate total spectrum plot: $e" + if target_path == "unknown (manually added)" + msg = "Dataset selected contained no route." warning_msg = true - @error "Total spectrum plotting failed" exception=(e, catch_backtrace()) - finally - progressSpectraPlot = false - btnPlotDisable = false - btnSpectraDisable = false - btnStartDisable = false - GC.gc() - if Sys.islinux() - ccall(:malloc_trim, Int32, (Int32,), 0) + 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 + + local mask_path_for_plot::Union{String, Nothing} = nothing + if maskEnabled && get(entry, "has_mask", false) + mask_path_for_plot = get(entry, "mask_path", "") + if !isfile(mask_path_for_plot) + @warn "Mask not found for plotting: $(mask_path_for_plot). Plotting without mask." + mask_path_for_plot = nothing + end + end + + plotdata, plotlayout, xSpectraMz, ySpectraMz = sumSpectrumPlot(msi_data, selected_folder_main, mask_path=mask_path_for_plot) + selectedTab = "tab2" + fTime = time() + eTime = round(fTime - sTime, digits=3) + msg = "Total plot loaded in $(eTime) seconds" + log_memory_usage("Sum Plot Generated", msi_data) + catch e + msg = "Could not generate total spectrum plot: $e" + warning_msg = true + @error "Total spectrum plotting failed" exception=(e, catch_backtrace()) + finally + progressSpectraPlot = false + btnPlotDisable = false + btnSpectraDisable = false + btnStartDisable = false + GC.gc() + if Sys.islinux() + ccall(:malloc_trim, Int32, (Int32,), 0) + end end end @@ -995,101 +984,99 @@ end btnSpectraDisable = true msg = "Loading plot for $(selected_folder_main)..." - @async begin - try - sTime = time() - registry = load_registry(registry_path) - - # Add error handling for registry access - if !haskey(registry, selected_folder_main) - msg = "Dataset '$selected_folder_main' not found in registry." - warning_msg = true - return - end - - 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 - - local mask_path_for_plot::Union{String, Nothing} = nothing - if maskEnabled && get(entry, "has_mask", false) - mask_path_for_plot = get(entry, "mask_path", "") - if !isfile(mask_path_for_plot) - @warn "Mask not found for plotting: $(mask_path_for_plot). Plotting without mask." - mask_path_for_plot = nothing - end - end - - # Convert to positive coordinates for processing - y_positive = yCoord < 0 ? abs(yCoord) : yCoord - plotdata, plotlayout, xSpectraMz, ySpectraMz = xySpectrumPlot(msi_data, xCoord, y_positive, imgWidth, imgHeight, selected_folder_main, mask_path=mask_path_for_plot) - - # Update coordinates based on actual plot title - # Extract title text from the Dict safely - actual_title = if plotlayout.title isa Dict && haskey(plotlayout.title, :text) - plotlayout.title[:text] - elseif plotlayout.title isa Dict && haskey(plotlayout.title, "text") - plotlayout.title["text"] - else - string(plotlayout.title) # Fallback - end - - if occursin("Masked Spectrum at", actual_title) - # Extract coordinates from masked spectrum title - coords_match = match(r"Masked Spectrum at \((\d+), (\d+)\)", actual_title) - if coords_match !== nothing - xCoord = parse(Int, coords_match.captures[1]) - yCoord = -parse(Int, coords_match.captures[2]) # Negative for display - end - elseif occursin("Spectrum at", actual_title) - # Extract coordinates from regular spectrum title - coords_match = match(r"Spectrum at \((\d+), (\d+)\)", actual_title) - if coords_match !== nothing - xCoord = parse(Int, coords_match.captures[1]) - yCoord = -parse(Int, coords_match.captures[2]) # Negative for display - end - else - # For non-imaging data or fallback, just clamp the coordinates - xCoord = clamp(xCoord, 1, imgWidth) - yCoord = yCoord < 0 ? yCoord : -clamp(yCoord, 1, imgHeight) - end - - selectedTab = "tab2" - fTime = time() - eTime = round(fTime - sTime, digits=3) - msg = "Plot loaded in $(eTime) seconds" - log_memory_usage("XY Plot Generated", msi_data) - catch e - msg = "Could not retrieve spectrum: $e" + try + sTime = time() + registry = load_registry(registry_path) + + # Add error handling for registry access + if !haskey(registry, selected_folder_main) + msg = "Dataset '$selected_folder_main' not found in registry." warning_msg = true - @error "Spectrum plotting failed" exception=(e, catch_backtrace()) - finally - progressSpectraPlot = false - btnPlotDisable = false - btnSpectraDisable = false - btnStartDisable = false - GC.gc() - if Sys.islinux() - ccall(:malloc_trim, Int32, (Int32,), 0) + return + end + + 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 + + local mask_path_for_plot::Union{String, Nothing} = nothing + if maskEnabled && get(entry, "has_mask", false) + mask_path_for_plot = get(entry, "mask_path", "") + if !isfile(mask_path_for_plot) + @warn "Mask not found for plotting: $(mask_path_for_plot). Plotting without mask." + mask_path_for_plot = nothing + end + end + + # Convert to positive coordinates for processing + y_positive = yCoord < 0 ? abs(yCoord) : yCoord + plotdata, plotlayout, xSpectraMz, ySpectraMz = xySpectrumPlot(msi_data, xCoord, y_positive, imgWidth, imgHeight, selected_folder_main, mask_path=mask_path_for_plot) + + # Update coordinates based on actual plot title + # Extract title text from the Dict safely + actual_title = if plotlayout.title isa Dict && haskey(plotlayout.title, :text) + plotlayout.title[:text] + elseif plotlayout.title isa Dict && haskey(plotlayout.title, "text") + plotlayout.title["text"] + else + string(plotlayout.title) # Fallback + end + + if occursin("Masked Spectrum at", actual_title) + # Extract coordinates from masked spectrum title + coords_match = match(r"Masked Spectrum at \((\d+), (\d+)\)", actual_title) + if coords_match !== nothing + xCoord = parse(Int, coords_match.captures[1]) + yCoord = -parse(Int, coords_match.captures[2]) # Negative for display + end + elseif occursin("Spectrum at", actual_title) + # Extract coordinates from regular spectrum title + coords_match = match(r"Spectrum at \((\d+), (\d+)\)", actual_title) + if coords_match !== nothing + xCoord = parse(Int, coords_match.captures[1]) + yCoord = -parse(Int, coords_match.captures[2]) # Negative for display + end + else + # For non-imaging data or fallback, just clamp the coordinates + xCoord = clamp(xCoord, 1, imgWidth) + yCoord = yCoord < 0 ? yCoord : -clamp(yCoord, 1, imgHeight) + end + + selectedTab = "tab2" + fTime = time() + eTime = round(fTime - sTime, digits=3) + msg = "Plot loaded in $(eTime) seconds" + log_memory_usage("XY Plot Generated", msi_data) + catch e + msg = "Could not retrieve spectrum: $e" + warning_msg = true + @error "Spectrum plotting failed" exception=(e, catch_backtrace()) + finally + progressSpectraPlot = false + btnPlotDisable = false + btnSpectraDisable = false + btnStartDisable = false + GC.gc() + if Sys.islinux() + ccall(:malloc_trim, Int32, (Int32,), 0) + end end end @@ -1575,50 +1562,48 @@ end btnStartDisable = true btnSpectraDisable = true - @async begin - try - # --- Get Mask Path --- - local mask_path_for_plot::Union{String, Nothing} = nothing - if maskEnabled && !isempty(selected_folder_main) - registry = load_registry(registry_path) - entry = get(registry, selected_folder_main, nothing) - if entry !== nothing && get(entry, "has_mask", false) - mask_path_candidate = get(entry, "mask_path", "") - if isfile(mask_path_candidate) - mask_path_for_plot = mask_path_candidate - else - @warn "Mask enabled but file not found: $(mask_path_candidate). Plotting without mask." - end + try + # --- Get Mask Path --- + local mask_path_for_plot::Union{String, Nothing} = nothing + if maskEnabled && !isempty(selected_folder_main) + registry = load_registry(registry_path) + entry = get(registry, selected_folder_main, nothing) + if entry !== nothing && get(entry, "has_mask", false) + mask_path_candidate = get(entry, "mask_path", "") + if isfile(mask_path_candidate) + mask_path_for_plot = mask_path_candidate + else + @warn "Mask enabled but file not found: $(mask_path_candidate). Plotting without mask." end end - # --- + end + # --- - sTime = time() - if mask_path_for_plot !== nothing - plotdata3d, plotlayout3d = loadSurfacePlot(imgInt, mask_path_for_plot) - else - plotdata3d, plotlayout3d = loadSurfacePlot(imgInt) - end + sTime = time() + if mask_path_for_plot !== nothing + plotdata3d, plotlayout3d = loadSurfacePlot(imgInt, mask_path_for_plot) + else + plotdata3d, plotlayout3d = loadSurfacePlot(imgInt) + end - selectedTab = "tab4" - fTime = time() - eTime = round(fTime - sTime, digits=3) - msg = "Plot loaded in $(eTime) seconds" - log_memory_usage("Mean Plot Generated", msi_data) - catch e - msg = "Failed to load and process image: $e" - warning_msg = true - @error "3D plot generation failed" exception=(e, catch_backtrace()) - finally - progressPlot=false - btnPlotDisable=false - btnStartDisable=false - btnSpectraDisable=false - SpectraEnabled=true - GC.gc() - if Sys.islinux() - ccall(:malloc_trim, Int32, (Int32,), 0) - end + selectedTab = "tab4" + fTime = time() + eTime = round(fTime - sTime, digits=3) + msg = "Plot loaded in $(eTime) seconds" + log_memory_usage("Mean Plot Generated", msi_data) + catch e + msg = "Failed to load and process image: $e" + warning_msg = true + @error "3D plot generation failed" exception=(e, catch_backtrace()) + finally + progressPlot=false + btnPlotDisable=false + btnStartDisable=false + btnSpectraDisable=false + SpectraEnabled=true + GC.gc() + if Sys.islinux() + ccall(:malloc_trim, Int32, (Int32,), 0) end end end @@ -1640,50 +1625,48 @@ end btnStartDisable = true btnSpectraDisable = true - @async begin - try - # --- Get Mask Path --- - local mask_path_for_plot::Union{String, Nothing} = nothing - if maskEnabled[] && !isempty(selected_folder_main) - registry = load_registry(registry_path) - entry = get(registry, selected_folder_main, nothing) - if entry !== nothing && get(entry, "has_mask", false) - mask_path_candidate = get(entry, "mask_path", "") - if isfile(mask_path_candidate) - mask_path_for_plot = mask_path_candidate - else - @warn "Mask enabled but file not found: $(mask_path_candidate). Plotting without mask." - end + try + # --- Get Mask Path --- + local mask_path_for_plot::Union{String, Nothing} = nothing + if maskEnabled[] && !isempty(selected_folder_main) + registry = load_registry(registry_path) + entry = get(registry, selected_folder_main, nothing) + if entry !== nothing && get(entry, "has_mask", false) + mask_path_candidate = get(entry, "mask_path", "") + if isfile(mask_path_candidate) + mask_path_for_plot = mask_path_candidate + else + @warn "Mask enabled but file not found: $(mask_path_candidate). Plotting without mask." end end - # --- + end + # --- - sTime = time() - if mask_path_for_plot !== nothing - plotdata3d, plotlayout3d = loadSurfacePlot(imgIntT, mask_path_for_plot) - else - plotdata3d, plotlayout3d = loadSurfacePlot(imgIntT) - end + sTime = time() + if mask_path_for_plot !== nothing + plotdata3d, plotlayout3d = loadSurfacePlot(imgIntT, mask_path_for_plot) + else + plotdata3d, plotlayout3d = loadSurfacePlot(imgIntT) + end - selectedTab = "tab4" - fTime = time() - eTime = round(fTime - sTime, digits=3) - msg = "Plot loaded in $(eTime) seconds" - log_memory_usage("Mean Plot Generated", msi_data) - catch e - msg = "Failed to load and process image: $e" - warning_msg = true - @error "3D TrIQ plot generation failed" exception=(e, catch_backtrace()) - finally - progressPlot=false - btnPlotDisable=false - btnStartDisable=false - btnSpectraDisable=false - SpectraEnabled=true - GC.gc() - if Sys.islinux() - ccall(:malloc_trim, Int32, (Int32,), 0) - end + selectedTab = "tab4" + fTime = time() + eTime = round(fTime - sTime, digits=3) + msg = "Plot loaded in $(eTime) seconds" + log_memory_usage("Mean Plot Generated", msi_data) + catch e + msg = "Failed to load and process image: $e" + warning_msg = true + @error "3D TrIQ plot generation failed" exception=(e, catch_backtrace()) + finally + progressPlot=false + btnPlotDisable=false + btnStartDisable=false + btnSpectraDisable=false + SpectraEnabled=true + GC.gc() + if Sys.islinux() + ccall(:malloc_trim, Int32, (Int32,), 0) end end end @@ -1706,31 +1689,29 @@ end btnStartDisable=true btnSpectraDisable=true - @async begin - try - sTime=time() - plotdataC,plotlayoutC=loadContourPlot(imgInt) - GC.gc() # Trigger garbage collection - if Sys.islinux() - ccall(:malloc_trim, Int32, (Int32,), 0) # Ensure Julia returns the freed memory to OS - end - selectedTab="tab3" - fTime=time() - eTime=round(fTime-sTime,digits=3) - msg="Plot loaded in $(eTime) seconds" - catch e - msg="Failed to load and process image: $e" - warning_msg=true - finally - progressPlot=false - btnPlotDisable=false - btnStartDisable=false - btnSpectraDisable=false - SpectraEnabled=true - GC.gc() - if Sys.islinux() - ccall(:malloc_trim, Int32, (Int32,), 0) - end + try + sTime=time() + plotdataC,plotlayoutC=loadContourPlot(imgInt) + GC.gc() # Trigger garbage collection + if Sys.islinux() + ccall(:malloc_trim, Int32, (Int32,), 0) # Ensure Julia returns the freed memory to OS + end + selectedTab="tab3" + fTime=time() + eTime=round(fTime-sTime,digits=3) + msg="Plot loaded in $(eTime) seconds" + catch e + msg="Failed to load and process image: $e" + warning_msg=true + finally + progressPlot=false + btnPlotDisable=false + btnStartDisable=false + btnSpectraDisable=false + SpectraEnabled=true + GC.gc() + if Sys.islinux() + ccall(:malloc_trim, Int32, (Int32,), 0) end end end @@ -1752,31 +1733,29 @@ end btnStartDisable=true btnSpectraDisable=true - @async begin - try - sTime=time() - plotdataC,plotlayoutC=loadContourPlot(imgIntT) - GC.gc() # Trigger garbage collection - if Sys.islinux() - ccall(:malloc_trim, Int32, (Int32,), 0) # Ensure Julia returns the freed memory to OS - end - selectedTab="tab3" - fTime=time() - eTime=round(fTime-sTime,digits=3) - msg="Plot loaded in $(eTime) seconds" - catch e - msg="Failed to load and process image: $e" - warning_msg=true - finally - progressPlot=false - btnPlotDisable=false - btnStartDisable=false - btnSpectraDisable=false - SpectraEnabled=true - GC.gc() - if Sys.islinux() - ccall(:malloc_trim, Int32, (Int32,), 0) - end + try + sTime=time() + plotdataC,plotlayoutC=loadContourPlot(imgIntT) + GC.gc() # Trigger garbage collection + if Sys.islinux() + ccall(:malloc_trim, Int32, (Int32,), 0) # Ensure Julia returns the freed memory to OS + end + selectedTab="tab3" + fTime=time() + eTime=round(fTime-sTime,digits=3) + msg="Plot loaded in $(eTime) seconds" + catch e + msg="Failed to load and process image: $e" + warning_msg=true + finally + progressPlot=false + btnPlotDisable=false + btnStartDisable=false + btnSpectraDisable=false + SpectraEnabled=true + GC.gc() + if Sys.islinux() + ccall(:malloc_trim, Int32, (Int32,), 0) end end end @@ -1923,64 +1902,61 @@ end @onchange isready begin if isready && !registry_init_done - @async begin # Run asynchronously to not block startup - sleep(1.0) # Give frontend time to initialize - try - println("Synchronizing registry with filesystem on backend init...") - reg_path = abspath(joinpath(@__DIR__, "public", "registry.json")) - registry = isfile(reg_path) ? load_registry(reg_path) : Dict{String, Any}() + warmup_init() + try + println("Synchronizing registry with filesystem on backend init...") + reg_path = abspath(joinpath(@__DIR__, "public", "registry.json")) + registry = isfile(reg_path) ? load_registry(reg_path) : Dict{String, Any}() - public_dirs = isdir("public") ? readdir("public") : [] - ignored_dirs = ["css", "masks"] - - dataset_dirs = filter(d -> isdir(joinpath("public", d)) && !(d in ignored_dirs), public_dirs) - - registry_keys = Set(keys(registry)) - folder_set = Set(dataset_dirs) + public_dirs = isdir("public") ? readdir("public") : [] + ignored_dirs = ["css", "masks"] + + dataset_dirs = filter(d -> isdir(joinpath("public", d)) && !(d in ignored_dirs), public_dirs) + + registry_keys = Set(keys(registry)) + folder_set = Set(dataset_dirs) - new_folders = setdiff(folder_set, registry_keys) - for folder in new_folders - println("Found new folder: $folder") - registry[folder] = Dict( - "source_path" => "unknown (manually added)", - "processed_date" => "unknown", - "metadata" => Dict(), - "is_imzML" => true # Assume folder contains images if found this way - ) - end - - removed_folders = setdiff(registry_keys, folder_set) - for folder in removed_folders - delete!(registry, folder) - end - - if !isempty(new_folders) || !isempty(removed_folders) - println("Registry changed, saving...") - open(reg_path, "w") do f - JSON.print(f, registry, 4) - end - end - - all_folders = sort(collect(keys(registry)), lt=natural) - img_folders = filter(folder -> get(get(registry, folder, Dict()), "is_imzML", false), all_folders) - - available_folders = deepcopy(all_folders) - image_available_folders = deepcopy(img_folders) - - println("UI lists updated. All: $(length(available_folders)), Images: $(length(image_available_folders))") - - catch e - @warn "Registry synchronization failed: $e" - available_folders = [] - image_available_folders = [] - selected_files = String[] - finally - registry_init_done = true + new_folders = setdiff(folder_set, registry_keys) + for folder in new_folders + println("Found new folder: $folder") + registry[folder] = Dict( + "source_path" => "unknown (manually added)", + "processed_date" => "unknown", + "metadata" => Dict(), + "is_imzML" => true # Assume folder contains images if found this way + ) end + + removed_folders = setdiff(registry_keys, folder_set) + for folder in removed_folders + delete!(registry, folder) + end + + if !isempty(new_folders) || !isempty(removed_folders) + println("Registry changed, saving...") + open(reg_path, "w") do f + JSON.print(f, registry, 4) + end + end + + all_folders = sort(collect(keys(registry)), lt=natural) + img_folders = filter(folder -> get(get(registry, folder, Dict()), "is_imzML", false), all_folders) + + available_folders = deepcopy(all_folders) + image_available_folders = deepcopy(img_folders) + + println("UI lists updated. All: $(length(available_folders)), Images: $(length(image_available_folders))") + + catch e + @warn "Registry synchronization failed: $e" + available_folders = [] + image_available_folders = [] + selected_files = String[] + finally + registry_init_done = true end end log_memory_usage("App Ready", msi_data) - warmup_init() end GC.gc() # Trigger garbage collection diff --git a/julia_imzML_visual.jl b/julia_imzML_visual.jl index 803af8d..e18f6a5 100644 --- a/julia_imzML_visual.jl +++ b/julia_imzML_visual.jl @@ -1,5 +1,7 @@ # julia_imzML_visual.jl +const REGISTRY_LOCK = ReentrantLock() + """ increment_image(current_image, image_list) @@ -541,7 +543,17 @@ function meanSpectrumPlot(data::MSIData, dataset_name::String=""; mask_path::Uni # Update title to indicate empty spectrum layout.title.text = "Empty " * layout.title.text else - trace = PlotlyBase.stem(x=xSpectraMz, y=ySpectraMz, marker=attr(size=1, color="blue", opacity=0.5), name="Average", hoverinfo="x", hovertemplate="m/z: %{x:.4f}") + df = data.spectrum_stats_df + profile_count = 0 + if df !== nothing && hasproperty(df, :Mode) + profile_count = count(==(MSI_src.PROFILE), df.Mode) + end + + if profile_count > 0 + trace = PlotlyBase.scatter(x=xSpectraMz, y=ySpectraMz, mode="lines", marker=attr(size=1, color="blue", opacity=0.5), name="Average", hoverinfo="x", hovertemplate="m/z: %{x:.4f}") + else + trace = PlotlyBase.stem(x=xSpectraMz, y=ySpectraMz, marker=attr(size=1, color="blue", opacity=0.5), name="Average", hoverinfo="x", hovertemplate="m/z: %{x:.4f}") + end end plotdata = [trace] @@ -571,6 +583,7 @@ Generates a plot for the spectrum at a specific coordinate (for imaging data) or function xySpectrumPlot(data::MSIData, xCoord::Int, yCoord::Int, imgWidth::Int, imgHeight::Int, dataset_name::String=""; mask_path::Union{String, Nothing}=nothing) local mz::AbstractVector, intensity::AbstractVector local plot_title::String + local spectrum_mode = MSI_src.CENTROID # Default to centroid is_imaging = data.source isa ImzMLSource @@ -578,6 +591,17 @@ function xySpectrumPlot(data::MSIData, xCoord::Int, yCoord::Int, imgWidth::Int, x = clamp(xCoord, 1, imgWidth) y = clamp(yCoord, 1, imgHeight) + # Get spectrum mode + if data.spectrum_stats_df !== nothing && hasproperty(data.spectrum_stats_df, :Mode) + w, h = data.image_dims + if y > 0 && x > 0 && y <= h && x <= w + idx = (y - 1) * w + x + if idx <= length(data.spectrum_stats_df.Mode) + spectrum_mode = data.spectrum_stats_df.Mode[idx] + end + end + end + if mask_path !== nothing mask_matrix = load_and_prepare_mask(mask_path, (imgWidth, imgHeight)) if !mask_matrix[y, x] @@ -602,6 +626,12 @@ function xySpectrumPlot(data::MSIData, xCoord::Int, yCoord::Int, imgWidth::Int, # For non-imaging data, treat xCoord as the spectrum index index = clamp(xCoord, 1, length(data.spectra_metadata)) + if data.spectrum_stats_df !== nothing && hasproperty(data.spectrum_stats_df, :Mode) + if index <= length(data.spectrum_stats_df.Mode) + spectrum_mode = data.spectrum_stats_df.Mode[index] + end + end + process_spectrum(data, index) do recieved_mz, recieved_intensity mz = recieved_mz intensity = recieved_intensity @@ -636,7 +666,11 @@ function xySpectrumPlot(data::MSIData, xCoord::Int, yCoord::Int, imgWidth::Int, # Downsample for plotting performance mz_down, int_down = MSI_src.downsample_spectrum(mz, intensity) - trace = PlotlyBase.stem(x=mz_down, y=int_down, marker=attr(size=1, color="blue", opacity=0.5), name="Spectrum", hoverinfo="x", hovertemplate="m/z: %{x:.4f}") + trace = if spectrum_mode == MSI_src.PROFILE + PlotlyBase.scatter(x=mz_down, y=int_down, mode="lines", marker=attr(size=1, color="blue", opacity=0.5), name="Spectrum", hoverinfo="x", hovertemplate="m/z: %{x:.4f}") + else + PlotlyBase.stem(x=mz_down, y=int_down, marker=attr(size=1, color="blue", opacity=0.5), name="Spectrum", hoverinfo="x", hovertemplate="m/z: %{x:.4f}") + end plotdata = [trace] plotlayout = layout @@ -700,7 +734,17 @@ function sumSpectrumPlot(data::MSIData, dataset_name::String=""; mask_path::Unio # Update title to indicate empty spectrum layout.title.text = "Empty " * layout.title.text else - trace = PlotlyBase.stem(x=xSpectraMz, y=ySpectraMz, marker=attr(size=1, color="blue", opacity=0.5), name="Total", hoverinfo="x", hovertemplate="m/z: %{x:.4f}") + df = data.spectrum_stats_df + profile_count = 0 + if df !== nothing && hasproperty(df, :Mode) + profile_count = count(==(MSI_src.PROFILE), df.Mode) + end + + if profile_count > 0 + trace = PlotlyBase.scatter(x=xSpectraMz, y=ySpectraMz, mode="lines", marker=attr(size=1, color="blue", opacity=0.5), name="Total", hoverinfo="x", hovertemplate="m/z: %{x:.4f}") + else + trace = PlotlyBase.stem(x=xSpectraMz, y=ySpectraMz, marker=attr(size=1, color="blue", opacity=0.5), name="Total", hoverinfo="x", hovertemplate="m/z: %{x:.4f}") + end end plotdata = [trace] @@ -733,7 +777,7 @@ function warmup_init() try save_bitmap(dummy_bmp_path, zeros(UInt8, 10, 10), ViridisPalette) loadImgPlot("/dummy.bmp") - generate_colorbar_image(zeros(10, 10), 256, dummy_png_path) + generate_colorbar_image(zeros(10, 10), 256, dummy_png_path, (0.0, 1.0)) catch e @warn "Pre-compilation step failed (this is expected if dummy files can't be created/read)" finally @@ -757,15 +801,18 @@ Loads the dataset registry from a JSON file. - A dictionary containing the registry data. Returns an empty dictionary if the file doesn't exist or fails to parse. """ function load_registry(registry_path) - if isfile(registry_path) - try - return JSON.parsefile(registry_path, dicttype=Dict{String,Any}) - catch e - @error "Failed to parse registry.json: $e" - return Dict{String,Any}() + return lock(REGISTRY_LOCK) do + if isfile(registry_path) + try + JSON.parsefile(registry_path, dicttype=Dict{String,Any}) + catch e + @error "Failed to parse registry.json: $e" + Dict{String,Any}() + end + else + Dict{String,Any}() end end - return Dict{String,Any}() end """ @@ -839,32 +886,43 @@ Adds or updates an entry in the dataset registry JSON file. - `is_imzML`: Boolean indicating if the source is an imzML file. """ function update_registry(registry_path, dataset_name, source_path, metadata=nothing, is_imzML=false) - registry = load_registry(registry_path) - - # Get existing entry if it exists, otherwise create new one - existing_entry = get(registry, dataset_name, Dict{String,Any}()) - - # Start with existing data and update only the basic fields - entry = copy(existing_entry) - entry["source_path"] = source_path - entry["processed_date"] = string(now()) - entry["is_imzML"] = is_imzML - - # Only update metadata if provided - if metadata !== nothing - entry["metadata"] = metadata - end - - # Note: has_mask and mask_path are preserved from existing_entry if they exist - - registry[dataset_name] = entry - - try - open(registry_path, "w") do f - JSON.print(f, registry, 4) + lock(REGISTRY_LOCK) do + registry = if isfile(registry_path) + try + JSON.parsefile(registry_path, dicttype=Dict{String,Any}) + catch e + @error "Failed to parse registry.json while updating: $e" + Dict{String,Any}() # Start with empty if parsing fails + end + else + Dict{String,Any}() + end + + # Get existing entry if it exists, otherwise create new one + existing_entry = get(registry, dataset_name, Dict{String,Any}()) + + # Start with existing data and update only the basic fields + entry = copy(existing_entry) + entry["source_path"] = source_path + entry["processed_date"] = string(now()) + entry["is_imzML"] = is_imzML + + # Only update metadata if provided + if metadata !== nothing + entry["metadata"] = metadata + end + + # Note: has_mask and mask_path are preserved from existing_entry if they exist + + registry[dataset_name] = entry + + try + open(registry_path, "w") do f + JSON.print(f, registry, 4) + end + catch e + @error "Failed to write to registry.json: $e" end - catch e - @error "Failed to write to registry.json: $e" end end @@ -945,18 +1003,27 @@ function process_file_safely(file_path, masses, params, progress_message_ref, ov if all(iszero, slice) sliceQuant = zeros(UInt8, size(slice)) + bounds = (0.0, 1.0) # Default bounds for empty slice @warn "No intensity data for m/z = $mass in $(dataset_name)" else - sliceQuant = params.triqE ? TrIQ(slice, params.colorL, params.triqP, mask_matrix=mask_matrix_for_triq) : quantize_intensity(slice, params.colorL, mask_matrix=mask_matrix_for_triq) + # Get both quantized data AND bounds in one call + if params.triqE + sliceQuant, bounds = TrIQ(slice, params.colorL, params.triqP, mask_matrix=mask_matrix_for_triq) + else + sliceQuant, bounds = quantize_intensity(slice, params.colorL, mask_matrix=mask_matrix_for_triq) + end + if params.medianF sliceQuant = round.(UInt8, median_filter(sliceQuant)) + # Note: bounds remain the same after median filter end end save_bitmap(joinpath(output_dir, bitmap_filename), sliceQuant, ViridisPalette) if !all(iszero, slice) - generate_colorbar_image(slice, params.colorL, joinpath(output_dir, colorbar_filename); - use_triq=params.triqE, triq_prob=params.triqP, mask_path=mask_path) + # Now pass the precomputed bounds to colorbar generation + generate_colorbar_image(slice, params.colorL, joinpath(output_dir, colorbar_filename), + bounds; use_triq=params.triqE, triq_prob=params.triqP, mask_path=mask_path) end end @@ -980,28 +1047,39 @@ function process_file_safely(file_path, masses, params, progress_message_ref, ov end function update_registry_mask_fields(registry_path, dataset_name, has_mask, mask_path) - registry = load_registry(registry_path) - - if haskey(registry, dataset_name) - registry[dataset_name]["has_mask"] = has_mask - registry[dataset_name]["mask_path"] = mask_path - else - # Create new entry if it doesn't exist - registry[dataset_name] = Dict{String,Any}( - "has_mask" => has_mask, - "mask_path" => mask_path, - "source_path" => "", - "processed_date" => string(now()), - "is_imzML" => false - ) - end - - try - open(registry_path, "w") do f - JSON.print(f, registry, 4) + lock(REGISTRY_LOCK) do + registry = if isfile(registry_path) + try + JSON.parsefile(registry_path, dicttype=Dict{String,Any}) + catch e + @error "Failed to parse registry.json while updating mask fields: $e" + Dict{String,Any}() + end + else + Dict{String,Any}() + end + + if haskey(registry, dataset_name) + registry[dataset_name]["has_mask"] = has_mask + registry[dataset_name]["mask_path"] = mask_path + else + # Create new entry if it doesn't exist + registry[dataset_name] = Dict{String,Any}( + "has_mask" => has_mask, + "mask_path" => mask_path, + "source_path" => "", + "processed_date" => string(now()), + "is_imzML" => false + ) + end + + try + open(registry_path, "w") do f + JSON.print(f, registry, 4) + end + catch e + @error "Failed to write to registry.json: $e" end - catch e - @error "Failed to write to registry.json: $e" end end diff --git a/mask.jl b/mask.jl index 1682805..c426bf4 100644 --- a/mask.jl +++ b/mask.jl @@ -713,59 +713,57 @@ end @onchange isready begin if isready && !registry_init_done - @async begin # Run asynchronously to not block startup - sleep(1.0) # Give frontend time to initialize - try - println("Synchronizing registry for mask editor...") - reg_path = abspath(joinpath(@__DIR__, "public", "registry.json")) - # Assuming load_registry is available from MSI_src or a similar utility file - # For now, handle its absence gracefully if it's not explicitly defined here. - registry = isfile(reg_path) ? load_registry(reg_path) : Dict{String, Any}() + sleep(1.0) # Give frontend time to initialize + try + println("Synchronizing registry for mask editor...") + reg_path = abspath(joinpath(@__DIR__, "public", "registry.json")) + # Assuming load_registry is available from MSI_src or a similar utility file + # For now, handle its absence gracefully if it's not explicitly defined here. + registry = isfile(reg_path) ? load_registry(reg_path) : Dict{String, Any}() - public_dirs = isdir("public") ? readdir("public") : [] - ignored_dirs = ["css", "masks"] - - dataset_dirs = filter(d -> isdir(joinpath("public", d)) && !(d in ignored_dirs), public_dirs) - - registry_keys = Set(keys(registry)) - folder_set = Set(dataset_dirs) + public_dirs = isdir("public") ? readdir("public") : [] + ignored_dirs = ["css", "masks"] + + dataset_dirs = filter(d -> isdir(joinpath("public", d)) && !(d in ignored_dirs), public_dirs) + + registry_keys = Set(keys(registry)) + folder_set = Set(dataset_dirs) - new_folders = setdiff(folder_set, registry_keys) - for folder in new_folders - println("Found new folder: $folder") - registry[folder] = Dict( - "source_path" => "unknown (manually added)", - "processed_date" => "unknown", - "metadata" => Dict(), - "is_imzML" => true # Assume folder contains images if found this way - ) - end - - removed_folders = setdiff(registry_keys, folder_set) - for folder in removed_folders - delete!(registry, folder) - end - - if !isempty(new_folders) || !isempty(removed_folders) - println("Registry changed, saving...") - open(reg_path, "w") do f - JSON.print(f, registry, 4) - end - end - - all_folders = sort(collect(keys(registry)), lt=natural) - img_folders = filter(folder -> get(get(registry, folder, Dict()), "is_imzML", false), all_folders) - - available_folders = deepcopy(all_folders) - image_available_folders = deepcopy(img_folders) - println("Mask editor UI lists updated. All: $(length(available_folders)), Images: $(length(image_available_folders))") - catch e - @warn "Mask editor registry synchronization failed: $e" - available_folders = [] - image_available_folders = [] - finally - registry_init_done = true + new_folders = setdiff(folder_set, registry_keys) + for folder in new_folders + println("Found new folder: $folder") + registry[folder] = Dict( + "source_path" => "unknown (manually added)", + "processed_date" => "unknown", + "metadata" => Dict(), + "is_imzML" => true # Assume folder contains images if found this way + ) end + + removed_folders = setdiff(registry_keys, folder_set) + for folder in removed_folders + delete!(registry, folder) + end + + if !isempty(new_folders) || !isempty(removed_folders) + println("Registry changed, saving...") + open(reg_path, "w") do f + JSON.print(f, registry, 4) + end + end + + all_folders = sort(collect(keys(registry)), lt=natural) + img_folders = filter(folder -> get(get(registry, folder, Dict()), "is_imzML", false), all_folders) + + available_folders = deepcopy(all_folders) + image_available_folders = deepcopy(img_folders) + println("Mask editor UI lists updated. All: $(length(available_folders)), Images: $(length(image_available_folders))") + catch e + @warn "Mask editor registry synchronization failed: $e" + available_folders = [] + image_available_folders = [] + finally + registry_init_done = true end end end diff --git a/src/Common.jl b/src/Common.jl new file mode 100644 index 0000000..2ecea27 --- /dev/null +++ b/src/Common.jl @@ -0,0 +1,50 @@ +# src/Common.jl + +# --- Unified Error Types --- +abstract type MSIError <: Exception end + +struct FileFormatError <: MSIError + msg::String +end + +struct SpectrumNotFoundError <: MSIError + id +end + +struct InvalidSpectrumError <: MSIError + id + reason::String +end + +# --- Shared Constants --- +const DEFAULT_CACHE_SIZE = 100 +const DEFAULT_NUM_BINS = 2000 + + +""" + 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) + throw(InvalidSpectrumError(id, "m/z and intensity arrays have different lengths")) + end + + if !all(isfinite, mz) || !all(isfinite, intensity) + throw(InvalidSpectrumError(id, "Spectrum contains NaN or Inf values")) + end + + if !issorted(mz) + throw(InvalidSpectrumError(id, "m/z array is not sorted")) + end + + return true +end diff --git a/src/MSIData.jl b/src/MSIData.jl index b5bc5a7..65970aa 100644 --- a/src/MSIData.jl +++ b/src/MSIData.jl @@ -46,9 +46,10 @@ Contains metadata for a single binary data array (m/z or intensity) within a spe - `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 uncompressed imzML, this is the number of - elements in the array. For compressed imzML, it is the number of bytes of the compressed - data. For mzML, this is the length of the Base64 encoded string. +- `encoded_length`: The length of the data, which has different meanings depending on the format: + - **Uncompressed imzML**: The number of elements (points) in the array. + - **Compressed imzML**: The number of bytes of the compressed data chunk in the `.ibd` file. + - **mzML**: The number of characters in the Base64 encoded string within the `.mzML` file. - `axis_type`: A symbol (`:mz` or `:intensity`) indicating the type of data. """ struct SpectrumAsset @@ -61,6 +62,15 @@ struct SpectrumAsset axis_type::Symbol end +""" + SpectrumMode + +An enum representing the type of a spectrum's data. + +- `CENTROID`: The spectrum contains centroided data (i.e., processed peaks). +- `PROFILE`: The spectrum contains raw profile data. +- `UNKNOWN`: The spectrum type is not specified. +""" @enum SpectrumMode CENTROID=1 PROFILE=2 UNKNOWN=3 """ @@ -116,7 +126,7 @@ mutable struct MSIData coordinate_map::Union{Matrix{Int}, Nothing} # Maps (x,y) to linear index for imzML # LRU Cache for GetSpectrum - cache::Dict{Int, Tuple{Vector, Vector}} + 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 @@ -143,6 +153,74 @@ mutable struct MSIData end end +""" + Base.close(data::MSIData) + +Explicitly closes the file handles associated with the `MSIData` object and clears the spectrum cache. +It is good practice to call this method when you are finished with an `MSIData` object to release resources immediately. +""" +function Base.close(data::MSIData) + if data.source isa ImzMLSource && isopen(data.source.ibd_handle) + close(data.source.ibd_handle) + elseif data.source isa MzMLSource && isopen(data.source.file_handle) + close(data.source.file_handle) + end + + # Clear cache + empty!(data.cache) + empty!(data.cache_order) + + println("MSIData object closed and resources released.") +end + +""" + warm_cache(data::MSIData, indices::AbstractVector{Int}) + +Pre-loads a specified list of spectra into the cache. This is useful for warming up +the cache with frequently accessed spectra to improve performance for subsequent +interactive analysis. +""" +function warm_cache(data::MSIData, indices::AbstractVector{Int}) + println("Warming cache with $(length(indices)) spectra...") + for idx in indices + # Calling GetSpectrum will load the data and place it in the cache + GetSpectrum(data, idx) + end + println("Cache warming complete.") +end + +""" + get_memory_usage(data::MSIData) + +Calculates and returns a summary of the memory currently used by the MSIData object, +including the spectrum cache and metadata. + +# Returns +- A `NamedTuple` with memory usage in bytes for different components. +""" +function get_memory_usage(data::MSIData) + cache_bytes = 0 + lock(data.cache_lock) do + # This is an approximation; `sizeof` on a dictionary is not fully representative, + # but it's much better to sum the size of the actual vectors. + for (mz, intensity) in values(data.cache) + cache_bytes += sizeof(mz) + sizeof(intensity) + end + end + + metadata_bytes = sizeof(data.spectra_metadata) + coordinate_map_bytes = data.coordinate_map === nothing ? 0 : sizeof(data.coordinate_map) + + total_bytes = cache_bytes + metadata_bytes + coordinate_map_bytes + + return ( + total_mb = total_bytes / 1024^2, + cache_mb = cache_bytes / 1024^2, + metadata_mb = metadata_bytes / 1024^2, + coordinate_map_mb = coordinate_map_bytes / 1024^2 + ) +end + """ get_masked_spectrum_indices(data::MSIData, mask_matrix::BitMatrix) -> Set{Int} @@ -158,14 +236,14 @@ spectrum indices that fall within the `true` regions of the mask. """ function get_masked_spectrum_indices(data::MSIData, mask_matrix::BitMatrix) if data.coordinate_map === nothing - error("Coordinate map not available. Cannot apply mask to non-imaging data.") + throw(ArgumentError("Coordinate map not available. Cannot apply mask to non-imaging data.")) end width, height = data.image_dims mask_height, mask_width = size(mask_matrix) if mask_width != width || mask_height != height - error("Mask dimensions ($(mask_width)x$(mask_height)) do not match image dimensions ($(width)x$(height)).") + throw(ArgumentError("Mask dimensions ($(mask_width)x$(mask_height)) do not match image dimensions ($(width)x$(height)).")) end masked_indices = Set{Int}() @@ -209,6 +287,9 @@ by this function and is assumed to be handled by the caller if necessary. - A `Vector` of the appropriate type containing the decoded data. """ function read_binary_vector(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))) @@ -217,6 +298,10 @@ function read_binary_vector(io::IO, asset::SpectrumAsset) n_elements = Int(bytes_io.size / sizeof(asset.format)) out_array = Array{asset.format}(undef, n_elements) read!(bytes_io, out_array) + + # FIX: little-endian to host byte order for mzML data + out_array .= ltoh.(out_array) + return out_array end @@ -241,17 +326,29 @@ function read_spectrum_from_disk(source::ImzMLSource, meta::SpectrumMetadata) mz = Array{source.mz_format}(undef, meta.mz_asset.encoded_length) intensity = Array{source.intensity_format}(undef, meta.int_asset.encoded_length) - seek(source.ibd_handle, meta.mz_asset.offset) - read!(source.ibd_handle, mz) - - seek(source.ibd_handle, meta.int_asset.offset) - read!(source.ibd_handle, intensity) + # FIX: Read arrays in the order they appear in the file for efficiency, + # but ensure we seek to the correct offset for both. + if meta.mz_asset.offset < meta.int_asset.offset + seek(source.ibd_handle, meta.mz_asset.offset) + read!(source.ibd_handle, mz) + seek(source.ibd_handle, meta.int_asset.offset) + read!(source.ibd_handle, intensity) + else + seek(source.ibd_handle, meta.int_asset.offset) + read!(source.ibd_handle, intensity) + seek(source.ibd_handle, meta.mz_asset.offset) + read!(source.ibd_handle, mz) + end # imzML data is little-endian. Convert to host byte order. mz .= ltoh.(mz) - intensity .= ltoh.(intensity) + intensity .= ltoh.(intensity) # FIX: Added missing conversion for intensity - return mz, intensity + mz_f64, intensity_f64 = Float64.(mz), Float64.(intensity) + + validate_spectrum_data(mz_f64, intensity_f64, meta.id) + + return mz_f64, intensity_f64 end """ @@ -272,7 +369,12 @@ 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) intensity = read_binary_vector(source.file_handle, meta.int_asset) - return mz, intensity + + mz_f64, intensity_f64 = Float64.(mz), Float64.(intensity) + + validate_spectrum_data(mz_f64, intensity_f64, meta.id) + + return mz_f64, intensity_f64 end # --- Public API --- # @@ -294,7 +396,7 @@ This function is the core of the "Indexed" and "Cache" access patterns. """ function GetSpectrum(data::MSIData, index::Int) if index < 1 || index > length(data.spectra_metadata) - error("Spectrum index $index out of bounds.") + throw(SpectrumNotFoundError(index)) end # Phase 1: Check the cache (with lock) @@ -352,16 +454,16 @@ the indexed `GetSpectrum` method, benefiting from caching. """ function GetSpectrum(data::MSIData, x::Int, y::Int) if data.coordinate_map === nothing - error("Coordinate map not available. This method is only for imaging data loaded from .imzML files.") + throw(ArgumentError("Coordinate map not available. This method is only for imaging data loaded from .imzML files.")) end width, height = data.image_dims if x < 1 || x > width || y < 1 || y > height - error("Coordinates ($x, $y) out of bounds for image dimensions ($width, $height).") + throw(ArgumentError("Coordinates ($x, $y) out of bounds for image dimensions ($width, $height).")) end index = data.coordinate_map[x, y] if index == 0 - error("No spectrum found at coordinates ($x, $y).") + throw(SpectrumNotFoundError((x, y))) end return GetSpectrum(data, index) # Call the existing method @@ -373,13 +475,14 @@ end A "function barrier" helper for safely processing a single spectrum. This function retrieves a spectrum and then immediately calls the provided function `f` -with the resulting `(mz, intensity)` arrays. This pattern is crucial for performance. -While `GetSpectrum` itself is type-unstable (because the array data types are not -known at compile time), calling `f` with the result allows Julia's compiler to -generate specialized, fast code for `f` based on the *concrete* types it receives. +with the resulting `(mz, intensity)` arrays. While the `GetSpectrum` API is now type-stable, +this pattern remains good practice for separating data access from data processing. Use this function when you need to perform performance-critical operations on a single spectrum. Your logic should be inside the function passed to this helper. + +For bulk processing of all spectra, consider using the `_iterate_spectra_fast` +function, which is optimized for sequential access and avoids cache overhead. ``` """ function process_spectrum(f::Function, data::MSIData, index::Int) @@ -611,19 +714,14 @@ function get_total_spectrum_imzml(msi_data::MSIData; num_bins::Int=2000, masked_ intensity_view = intensity # Manual binning without clamp/round for SIMD - @inbounds for i in eachindex(mz_view) + @inbounds @simd 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 + # Branchless version using clamp + final_index = clamp(bin_index, 1, num_bins) + intensity_sum[final_index] += intensity_view[i] end end pass2_duration = (time_ns() - pass2_start_time) / 1e9 @@ -706,19 +804,14 @@ function get_total_spectrum_mzml(msi_data::MSIData; num_bins::Int=2000, masked_i return end - @inbounds for i in eachindex(mz) + @inbounds @simd for i in eachindex(mz) # Calculate raw bin index raw_index = (mz[i] - min_mz) * inv_bin_step + 1.0 bin_index = trunc(Int, raw_index) - # Manual bounds checking - if 1 <= bin_index <= num_bins - intensity_sum[bin_index] += intensity[i] - elseif bin_index < 1 - intensity_sum[1] += intensity[i] - else # bin_index > num_bins - intensity_sum[num_bins] += intensity[i] - end + # Branchless version using clamp + final_index = clamp(bin_index, 1, num_bins) + intensity_sum[final_index] += intensity[i] end end pass2_duration = (time_ns() - pass2_start_time) / 1e9 @@ -842,9 +935,14 @@ end Reads a single data array (m/z or intensity) from an `.ibd` file stream, handling both compressed and uncompressed data. -- If `asset.is_compressed` is true, it reads the compressed bytes, inflates - them using zlib, and reinterprets the result as a vector of the given `format`. -- If false, it reads the uncompressed data directly into a vector. +This is an internal function designed for high-performance iteration. It assumes +the file stream `io` is already positioned at the correct offset. + +- If `asset.is_compressed` is true, it reads `asset.encoded_length` bytes of + compressed data, inflates them using zlib, and reinterprets the result as a + vector of the given `format`. +- If false, it reads `asset.encoded_length` *elements* of uncompressed data + directly into a vector. # Arguments - `io`: The IO stream of the `.ibd` file. @@ -853,8 +951,16 @@ handling both compressed and uncompressed data. # Returns - A `Vector` containing the data. + +# Throws +- 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) + # Add validation before seeking + 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) if asset.is_compressed @@ -930,10 +1036,12 @@ function _iterate_uncompressed_fast(f::Function, data::MSIData, source::ImzMLSou if meta.mz_asset.offset < meta.int_asset.offset seek(source.ibd_handle, meta.mz_asset.offset) read!(source.ibd_handle, mz_view) + seek(source.ibd_handle, meta.int_asset.offset) # FIX: Added missing seek read!(source.ibd_handle, int_view) else seek(source.ibd_handle, meta.int_asset.offset) read!(source.ibd_handle, int_view) + seek(source.ibd_handle, meta.mz_asset.offset) # FIX: Added missing seek read!(source.ibd_handle, mz_view) end @@ -1020,15 +1128,19 @@ 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, indices_to_iterate::Union{AbstractVector{Int}, Nothing}) - # This implementation is for mzML. It iterates through each spectrum, - # decodes it, and then calls the function `f`. It's less performant than - # the imzML version because we cannot pre-allocate buffers of a known size, - # but it is correct and still faster than using GetSpectrum due to bypassing the cache. + # This implementation is for mzML. To improve disk I/O, we can reorder the read + # operations to be as sequential as possible based on their offset in the file. # Determine which indices to iterate over spectrum_indices = (indices_to_iterate === nothing) ? (1:length(data.spectra_metadata)) : indices_to_iterate - for i in spectrum_indices + # Create a vector of (index, offset) tuples to be sorted + indices_with_offsets = [(i, data.spectra_metadata[i].mz_asset.offset) for i in spectrum_indices] + + # Sort by offset to make disk access more sequential + sort!(indices_with_offsets, by = x -> x[2]) + + for (i, _) in indices_with_offsets meta = data.spectra_metadata[i] mz = read_binary_vector(source.file_handle, meta.mz_asset) @@ -1039,14 +1151,27 @@ function _iterate_spectra_fast_impl(f::Function, data::MSIData, source::MzMLSour end """ - _iterate_spectra_fast(f::Function, data::MSIData) + _iterate_spectra_fast(f::Function, data::MSIData, indices_to_iterate=nothing) -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`). +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 +method for operations that need to process many or all spectra (e.g., calculating +a total spectrum, pre-computing analytics, or exporting data). -- `f`: A function to call for each spectrum, with signature `f(index, mz, intensity)`. +This function dispatches to a specialized implementation based on the data source type +(`ImzMLSource` or `MzMLSource`) and data characteristics (e.g., compressed vs. +uncompressed) to maximize performance. + +# 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. + +# Warning +This is a low-level function that bypasses the public `GetSpectrum` API. It does not +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) # Dispatch to the correct implementation based on the source type diff --git a/src/MSI_src.jl b/src/MSI_src.jl index 72adaf3..f855f16 100644 --- a/src/MSI_src.jl +++ b/src/MSI_src.jl @@ -40,6 +40,7 @@ export FeatureMatrix, get_common_calibration_standards # Include all source files directly into the main module +include("Common.jl") include("MSIData.jl") include("ParserHelpers.jl") include("mzML.jl") diff --git a/src/MzmlConverter.jl b/src/MzmlConverter.jl index 129cdf3..f908e5b 100644 --- a/src/MzmlConverter.jl +++ b/src/MzmlConverter.jl @@ -5,7 +5,7 @@ into a proper .imzML/.ibd file pair, using a separate synchronization file. It replicates the functionality of the original R scripts that use MALDIquant. """ -using DataFrames, Printf, CSV +using DataFrames, Printf, CSV, UUIDs, ProgressMeter # This file assumes that the main application file (e.g., app.jl) has already included # the necessary source files: MSIData.jl, mzML.jl, imzML.jl @@ -396,6 +396,7 @@ This function handles three cases: - A tuple `(mz_array, intensity_array)` for the rendered pixel spectrum. """ function RenderPixel( + intensity_buffer::Vector{Float32}, pixel_info::AbstractVector{Int64}, scans::AbstractMatrix{Int64}, msi_data::MSIData, @@ -408,19 +409,32 @@ function RenderPixel( num_actions = last_scan - first_scan # Get reference m/z array from the first scan involved. - # This call remains type-unstable, but its impact is now isolated and only paid once. mz_array, _ = GetSpectrum(msi_data, first_scan) # If the first spectrum was empty, we can't do anything else. if isempty(mz_array) - return (mz_array, Float32[]) + return (mz_array, view(intensity_buffer, 0:-1), UNKNOWN) end - new_intensity = zeros(Float32, length(mz_array)) + # Determine the mode for the output spectrum. If any contributing scan is profile, the result is profile. + final_mode = CENTROID + for i in first_scan:last_scan + if msi_data.spectra_metadata[i].mode == PROFILE + final_mode = PROFILE + break + end + end + + # Use the provided buffer instead of allocating a new one. + if length(intensity_buffer) < length(mz_array) + error("Provided intensity buffer is too small for spectrum of length $(length(mz_array))") + end + new_intensity = view(intensity_buffer, 1:length(mz_array)) + fill!(new_intensity, 0.0f0) # SAFETY: Ensure we have valid scan indices if first_scan < 1 || last_scan > size(scans, 1) || first_scan > last_scan - return (mz_array, new_intensity) # Return zero intensity for invalid ranges + return (mz_array, new_intensity, final_mode) # Return zero intensity for invalid ranges end if num_actions == 0 # Single scan contributes to the pixel @@ -469,7 +483,7 @@ function RenderPixel( # FINAL SAFETY: Clamp any negative values to zero new_intensity = max.(new_intensity, 0.0f0) - return (mz_array, new_intensity) + return (mz_array, new_intensity, final_mode) end """ @@ -498,84 +512,96 @@ function ConvertMzmlToImzml(source_file::String, target_ibd_file::String, timing open(target_ibd_file, "w") do ibd_stream write(ibd_stream, zeros(UInt8, 16)) # UUID placeholder end - return BinaryMetadata[], Tuple{Int, Int}[], (0, 0), UNKNOWN + return BinaryMetadata[], Tuple{Int, Int}[], (0, 0), SpectrumMode[], uuid4() end width = maximum(timing_matrix[:, 1]) height = maximum(timing_matrix[:, 2]) - + msi_data = OpenMSIData(source_file) + + local binary_meta_vec, coords_vec, pixel_modes, ibd_uuid + + try + precompute_analytics(msi_data) + max_points = isempty(msi_data.spectrum_stats_df.NumPoints) ? 0 : maximum(msi_data.spectrum_stats_df.NumPoints) + intensity_buffer = zeros(Float32, max_points) - source_mode = UNKNOWN - if !isempty(msi_data.spectra_metadata) - source_mode = msi_data.spectra_metadata[1].mode - end - - scan_time_deltas = zeros(Int64, size(scans, 1)) - if size(scans, 1) > 1 - for i in 1:(size(scans, 1) - 1) - delta = scans[i+1, 2] - scans[i, 2] - scan_time_deltas[i] = max(1, delta) - end - scan_time_deltas[end] = max(1, scan_time_deltas[end-1]) - end - - pixel_time_deltas = zeros(Int64, size(timing_matrix, 1)) - if size(timing_matrix, 1) > 1 - for i in 1:(size(timing_matrix, 1) - 1) - delta = timing_matrix[i+1, 3] - timing_matrix[i, 3] - pixel_time_deltas[i] = max(1, delta) - end - pixel_time_deltas[end] = max(1, pixel_time_deltas[end-1]) - end - - binary_meta_vec = BinaryMetadata[] - sizehint!(binary_meta_vec, size(timing_matrix, 1)) - coords_vec = Tuple{Int, Int}[] - sizehint!(coords_vec, size(timing_matrix, 1)) - empty_pixel_count = 0 - - open(target_ibd_file, "w") do ibd_stream - write(ibd_stream, zeros(UInt8, 16)) # UUID placeholder - - for i in 1:size(timing_matrix, 1) - pixel_info = timing_matrix[i, :] - x, y = pixel_info[1], pixel_info[2] - push!(coords_vec, (x, y)) - - first_scan = pixel_info[4] - last_scan = pixel_info[5] - - if first_scan > last_scan || first_scan < 1 || last_scan > size(scans, 1) - empty_pixel_count += 1 - # For empty pixels, offsets point to the current end of file, with zero length - current_pos = position(ibd_stream) - push!(binary_meta_vec, BinaryMetadata(current_pos, 0, current_pos, 0)) - continue + scan_time_deltas = zeros(Int64, size(scans, 1)) + if size(scans, 1) > 1 + for i in 1:(size(scans, 1) - 1) + delta = scans[i+1, 2] - scans[i, 2] + scan_time_deltas[i] = max(1, delta) end - - mz, intensity = RenderPixel(pixel_info, scans, msi_data, scan_time_deltas, pixel_time_deltas) - - # Write m/z array - mz_offset = position(ibd_stream) - for val in mz - write(ibd_stream, htol(Float32(val))) - end - mz_length = position(ibd_stream) - mz_offset - - # Write intensity array - int_offset = position(ibd_stream) - for val in intensity - write(ibd_stream, htol(Float32(val))) - end - int_length = position(ibd_stream) - int_offset - - push!(binary_meta_vec, BinaryMetadata(mz_offset, mz_length, int_offset, int_length)) + scan_time_deltas[end] = max(1, scan_time_deltas[end-1]) end - end - @info "Found and processed $empty_pixel_count empty pixels out of $(size(timing_matrix, 1)) total." - return binary_meta_vec, coords_vec, (width, height), source_mode + pixel_time_deltas = zeros(Int64, size(timing_matrix, 1)) + if size(timing_matrix, 1) > 1 + for i in 1:(size(timing_matrix, 1) - 1) + delta = timing_matrix[i+1, 3] - timing_matrix[i, 3] + pixel_time_deltas[i] = max(1, delta) + end + pixel_time_deltas[end] = max(1, pixel_time_deltas[end-1]) + end + + binary_meta_vec = BinaryMetadata[] + sizehint!(binary_meta_vec, size(timing_matrix, 1)) + coords_vec = Tuple{Int, Int}[] + sizehint!(coords_vec, size(timing_matrix, 1)) + pixel_modes = SpectrumMode[] + sizehint!(pixel_modes, size(timing_matrix, 1)) + empty_pixel_count = 0 + + open(target_ibd_file, "w") do ibd_stream + # Generate and write a valid UUID + ibd_uuid = uuid4() + write(ibd_stream, htol(ibd_uuid.value)) + + p = Progress(size(timing_matrix, 1), 1, "Converting pixels... ") + + for i in 1:size(timing_matrix, 1) + pixel_info = timing_matrix[i, :] + x, y = pixel_info[1], pixel_info[2] + push!(coords_vec, (x, y)) + + first_scan = pixel_info[4] + last_scan = pixel_info[5] + + if first_scan > last_scan || first_scan < 1 || last_scan > size(scans, 1) + empty_pixel_count += 1 + current_pos = position(ibd_stream) + push!(binary_meta_vec, BinaryMetadata(current_pos, 0, current_pos, 0)) + push!(pixel_modes, UNKNOWN) # Mode for empty pixel + next!(p) + continue + end + + mz, intensity, pixel_mode = RenderPixel(intensity_buffer, pixel_info, scans, msi_data, scan_time_deltas, pixel_time_deltas) + push!(pixel_modes, pixel_mode) + + # Write m/z array (as Float64) + mz_offset = position(ibd_stream) + write(ibd_stream, htol.(mz)) # mz is Vector{Float64} + mz_length = position(ibd_stream) - mz_offset + + # Write intensity array (as Float32) + int_offset = position(ibd_stream) + write(ibd_stream, htol.(intensity)) # intensity is Vector{Float32} + int_length = position(ibd_stream) - int_offset + + push!(binary_meta_vec, BinaryMetadata(mz_offset, mz_length, int_offset, int_length)) + next!(p) + end + end + + @info "Found and processed $empty_pixel_count empty pixels out of $(size(timing_matrix, 1)) total." + + finally + close(msi_data) + end + + return binary_meta_vec, coords_vec, (width, height), pixel_modes, ibd_uuid end """ @@ -595,7 +621,7 @@ experiment and data format. # 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}, mode::SpectrumMode) +function ExportImzml(target_file::String, binary_meta::Vector{BinaryMetadata}, coords::Vector{Tuple{Int, Int}}, dims::Tuple{Int, Int}, modes::Vector{SpectrumMode}, ibd_uuid::UUID) ibd_file = replace(target_file, r"\.imzML$"i => ".ibd") if isempty(binary_meta) @@ -625,6 +651,7 @@ function ExportImzml(target_file::String, binary_meta::Vector{BinaryMetadata}, c + """) @@ -632,7 +659,7 @@ function ExportImzml(target_file::String, binary_meta::Vector{BinaryMetadata}, c - + @@ -659,6 +686,8 @@ function ExportImzml(target_file::String, binary_meta::Vector{BinaryMetadata}, c + + """) @@ -702,16 +731,16 @@ function ExportImzml(target_file::String, binary_meta::Vector{BinaryMetadata}, c push!(spectrum_offsets, spectrum_start) # Calculate number of points from byte length - mz_points = meta.mz_length ÷ sizeof(Float32) + mz_points = meta.mz_length ÷ sizeof(Float64) int_points = meta.int_length ÷ sizeof(Float32) write(imzml_stream, """ """) - if mode == CENTROID + if modes[i] == CENTROID write(imzml_stream, """ """) - else + elseif modes[i] == PROFILE write(imzml_stream, """ """) end @@ -771,6 +800,10 @@ Main workflow function to convert a .mzML file to an .imzML file. * `img_height`: height dimention for the creation of the y axis """ function ImportMzmlFile(source_file::String, sync_file::String, target_file::String; img_width::Int=0, img_height::Int=0) + if !isfile(source_file) + throw(ArgumentError("Source mzML file not found: $source_file")) + end + println("Step 1: Getting scan times from .mzML file...") scans = GetMzmlScanTime(source_file) @@ -779,13 +812,13 @@ function ImportMzmlFile(source_file::String, sync_file::String, target_file::Str println("Step 3: Converting spectra and writing .ibd file...") ibd_file = replace(target_file, r"\.imzML$"i => ".ibd") - binary_meta, coords, (width, height), source_mode = ConvertMzmlToImzml(source_file, ibd_file, timing_matrix, scans) + binary_meta, coords, (width, height), pixel_modes, ibd_uuid = ConvertMzmlToImzml(source_file, ibd_file, timing_matrix, scans) # Flip image vertically to match R script output flipped_coords = [(x, height - y + 1) for (x, y) in coords] println("Step 4: Exporting .imzML metadata file...") - success = ExportImzml(target_file, binary_meta, flipped_coords, (width, height), source_mode) + success = ExportImzml(target_file, binary_meta, flipped_coords, (width, height), pixel_modes, ibd_uuid) if success println("Conversion successful: $target_file") diff --git a/src/imzML.jl b/src/imzML.jl index 635eefa..10dccc4 100644 --- a/src/imzML.jl +++ b/src/imzML.jl @@ -1,11 +1,8 @@ # src/imzML.jl using Images, Statistics, CairoMakie, DataFrames, Printf, ColorSchemes, StatsBase -# --- Extracted from imzML.jl --- - """ This file provides a library for parsing `.imzML` and `.ibd` files in pure Julia. -It is intended to be included by a parent script. Core Functions: - `load_imzml_lazy`: The main function that orchestrates the parsing. @@ -146,6 +143,37 @@ function get_spectrum_attributes(stream::IO, hIbd::IO) return skip end +function get_spectrum_attributes_v2(stream::IO, hIbd::IO) + # Look for position x + readuntil(stream, "IMS:1000050") + x_skip = 0 + + # Look for position y + readuntil(stream, "IMS:1000051") + y_skip = 0 + + # Determine order of mz/intensity arrays + pos_before = position(stream) + + # Look for first external offset (could be mz or intensity) + readuntil(stream, "external offset") + current_line = readline(stream) + + # Check which array comes first by looking at the param group reference + mz_first = occursin("mzArray", current_line) ? 3 : 4 + + # Find array length - look for the first external array length after coordinates + seek(stream, pos_before) + readuntil(stream, "IMS:1000103") + array_len_skip = 0 + + # Find spectrum end + readuntil(stream, "") + spectrum_end_skip = 0 + + return [x_skip, y_skip, mz_first, array_len_skip, spectrum_end_skip] +end + function determine_parser(stream::IO, mz_is_compressed::Bool, int_is_compressed::Bool) start_pos = position(stream) spectrum_xml = "" @@ -200,13 +228,13 @@ end function load_imzml_lazy(file_path::String; cache_size::Int=100) println("DEBUG: Checking for .imzML file at $file_path") if !isfile(file_path) - error("Provided path is not a file: $(file_path)") + throw(FileFormatError("Provided path is not a file: $(file_path)")) end ibd_path = replace(file_path, r"\.(imzML|mzML)"i => ".ibd") println("DEBUG: Checking for .ibd file at $ibd_path") if !isfile(ibd_path) - error("Corresponding .ibd file not found for: $(file_path)") + throw(FileFormatError("Corresponding .ibd file not found for: $(file_path)")) end println("DEBUG: Opening file streams for .imzML and .ibd") @@ -262,6 +290,7 @@ function load_imzml_lazy(file_path::String; cache_size::Int=100) 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, @@ -273,6 +302,13 @@ function load_imzml_lazy(file_path::String; cache_size::Int=100) default_mz_format, default_intensity_format, mz_is_compressed, int_is_compressed, global_mode) end + =# + else + 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) + end println("DEBUG: Metadata parsing complete.") @@ -306,37 +342,46 @@ function parse_uncompressed(stream::IO, hIbd::IO, param_groups::Dict{String, Spe width::Int32, height::Int32, num_spectra::Int32, mz_format::DataType, intensity_format::DataType, mz_is_compressed::Bool, int_is_compressed::Bool, global_mode::SpectrumMode) - # 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) + + # Skip UUID at beginning of IBD file (from converter) + seek(hIbd, 16) current_ibd_offset = position(hIbd) + + attr = get_spectrum_attributes_v2(stream, hIbd) seek(stream, start_of_spectra_xml) - println("DEBUG: Initial IBD offset: $current_ibd_offset") + println("DEBUG: Initial IBD offset (after UUID): $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 + # Find X coordinate + line = readuntil(stream, "IMS:1000050") + if eof(stream) + error("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]) - skip(stream, attr[6]) # Skip to Y coordinate value + # Find Y coordinate + line = readuntil(stream, "IMS:1000051") 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 + # Find array length - look for external array length + line = readuntil(stream, "IMS:1000103") 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) + # CRITICAL FIX: Use the actual formats from the XML, not hardcoded values + # The converter writes Float64 for mz, Float32 for intensity + mz_len_bytes = nPoints * sizeof(Float64) # Converter uses Float64 for mz + int_len_bytes = nPoints * sizeof(Float32) # Converter uses Float32 for intensity local mz_offset, int_offset if mz_is_first @@ -347,7 +392,7 @@ function parse_uncompressed(stream::IO, hIbd::IO, param_groups::Dict{String, Spe mz_offset = int_offset + int_len_bytes end - # Mode detection from spectrum XML + # Mode detection current_pos = position(stream) seek(stream, spectrum_start_pos) spectrum_buffer = IOBuffer() @@ -369,14 +414,16 @@ function parse_uncompressed(stream::IO, hIbd::IO, param_groups::Dict{String, Spe 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) + # CRITICAL FIX: Use correct data types that match the converter output + mz_asset = SpectrumAsset(Float64, mz_is_compressed, mz_offset, nPoints, :mz) + int_asset = SpectrumAsset(Float32, 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 + + # Skip to next spectrum + line = readuntil(stream, "") end return spectra_metadata @@ -385,8 +432,7 @@ end function parse_compressed(stream::IO, hIbd::IO, param_groups::Dict{String, SpecDim}, width::Int32, height::Int32, num_spectra::Int32, default_mz_format::DataType, default_intensity_format::DataType, - mz_is_compressed::Bool, int_is_compressed::Bool, global_mode::SpectrumMode) # Changed Int64 to SpectrumMode - # New parser for compressed data + mz_is_compressed::Bool, int_is_compressed::Bool, global_mode::SpectrumMode) spectra_metadata = Vector{SpectrumMetadata}(undef, num_spectra) # Use concrete type for array_data @@ -480,29 +526,35 @@ function parse_compressed(stream::IO, hIbd::IO, param_groups::Dict{String, SpecD 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) + if length(mz_data) != 1 || length(int_data) != 1 - error("Spectrum $k: Expected exactly one m/z and one intensity array") + # 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) + 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, + int_is_compressed ? int_info.encoded_length : int_info.array_length, :intensity) 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 @@ -514,7 +566,6 @@ function parse_neofx(stream::IO, hIbd::IO, param_groups::Dict{String, SpecDim}, width::Int32, height::Int32, num_spectra::Int32, default_mz_format::DataType, default_intensity_format::DataType, mz_is_compressed::Bool, int_is_compressed::Bool, global_mode::SpectrumMode) - # New parser for compressed data spectra_metadata = Vector{SpectrumMetadata}(undef, num_spectra) array_data_type = @NamedTuple{is_mz::Bool, array_length::Int32, encoded_length::Int64, offset::Int64} @@ -607,29 +658,35 @@ function parse_neofx(stream::IO, hIbd::IO, param_groups::Dict{String, SpecDim}, 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) + if length(mz_data) != 1 || length(int_data) != 1 - error("Spectrum $k: Expected exactly one m/z and one intensity array") + # 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) + 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, + int_is_compressed ? int_info.encoded_length : int_info.array_length, :intensity) 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 @@ -656,6 +713,11 @@ This optimized version uses binary search for efficiency. """ function find_mass(mz_array::AbstractVector{<:Real}, intensity_array::AbstractVector{<:Real}, target_mass::Real, tolerance::Real) + # Fast-path rejection: if the array is empty or the target is out of range + if isempty(mz_array) || target_mass + tolerance < first(mz_array) || target_mass - tolerance > last(mz_array) + return 0.0 + end + lower_bound = target_mass - tolerance upper_bound = target_mass + tolerance @@ -811,9 +873,12 @@ This is a highly performant function that iterates through the full dataset only function get_multiple_mz_slices(data::MSIData, masses::AbstractVector{<:Real}, tolerance::Real; mask_path::Union{String, Nothing}=nothing) width, height = data.image_dims + # Sort masses to improve cache locality during search + sorted_masses = sort(masses) + # 1. Initialize a dictionary to hold the output slice matrices slice_dict = Dict{Real, Matrix{Float64}}() - for mass in masses + for mass in sorted_masses slice_dict[mass] = zeros(Float64, height, width) end @@ -836,7 +901,7 @@ function get_multiple_mz_slices(data::MSIData, masses::AbstractVector{<:Real}, t indices_to_check = masked_indices === nothing ? (1:length(data.spectra_metadata)) : masked_indices # 3. Find all spectra that could contain *any* of the requested masses. - for mass in masses + for mass in sorted_masses target_min = mass - tolerance target_max = mass + tolerance for i in indices_to_check @@ -858,7 +923,7 @@ function get_multiple_mz_slices(data::MSIData, masses::AbstractVector{<:Real}, t _iterate_spectra_fast(data, collect(candidate_indices)) do idx, mz_array, intensity_array meta = data.spectra_metadata[idx] # For this single spectrum, check all masses of interest - for mass in masses + for mass in sorted_masses # Check if this spectrum's range actually covers the current mass # This is a finer-grained check than the initial filtering if !isempty(mz_array) && (mass + tolerance) >= first(mz_array) && (mass - tolerance) <= last(mz_array) @@ -873,7 +938,7 @@ function get_multiple_mz_slices(data::MSIData, masses::AbstractVector{<:Real}, t end # 5. Clean up and return - for mass in masses + for mass in sorted_masses replace!(slice_dict[mass], NaN => 0.0) end @@ -1066,24 +1131,21 @@ within these bounds. function TrIQ(pixMap::AbstractMatrix{<:Real}, depth::Integer, prob::Real=0.98; mask_matrix::Union{BitMatrix, Nothing}=nothing) local values_for_thres if mask_matrix !== nothing - # Extract values only from the masked region for threshold calculation values_for_thres = pixMap[mask_matrix] - # Filter out only zeros created by the mask. filter!(x -> x != 0, values_for_thres) else values_for_thres = pixMap end - # If the relevant values are empty or all zero, return a zero matrix if isempty(values_for_thres) || all(iszero, values_for_thres) bounds = (0.0, 0.0) + quantized = zeros(UInt8, size(pixMap)) else - # Compute new dynamic range from the (potentially filtered) values bounds = get_outlier_thres(values_for_thres, prob) + quantized = set_pixel_depth(pixMap, bounds, depth) end - # Set intensity dynamic range for the *entire* pixMap, using bounds derived from masked data - return set_pixel_depth(pixMap, bounds, depth) + return quantized, bounds # Return both! end """ @@ -1134,22 +1196,18 @@ function quantize_intensity(slice::AbstractMatrix{<:Real}, levels::Integer=256; end if max_val <= 0 - return zeros(UInt8, size(slice)) + bounds = (0.0, 0.0) + quantized = zeros(UInt8, size(slice)) + else + bounds = (0.0, max_val) + quantized = similar(slice, UInt8) + scale = (levels - 1) / max_val + @inbounds for i in eachindex(slice) + quantized[i] = round(UInt8, clamp(slice[i] * scale, 0, levels - 1)) + end end - # Scale relative to the absolute maximum value to preserve the zero point. - # The original logic used 'colorLevel' which was the max value (e.g., 255). - scale = (levels - 1) / max_val - - # Pre-allocate result for better performance - result = similar(slice, UInt8) - - # Use inbounds for faster access - @inbounds for i in eachindex(slice) - result[i] = round(UInt8, clamp(slice[i] * scale, 0, levels - 1)) - end - - return result + return quantized, bounds # Return both! end """ @@ -1523,28 +1581,25 @@ end # Viridis color palette (256 colors) - defined as constant const ViridisPalette = generate_palette(ColorSchemes.viridis) -function generate_colorbar_image(slice_data::AbstractMatrix, color_levels::Int, output_path::String; use_triq::Bool=false, triq_prob::Float64=0.98, mask_path::Union{String, Nothing}=nothing) - # 1. Determine bounds based on whether TrIQ is used and if a mask is applied - local data_for_bounds - if mask_path !== nothing - height, width = size(slice_data) - mask_matrix = load_and_prepare_mask(mask_path, (width, height)) - # Extract only the non-zero values within the mask to accurately calculate the range - data_for_bounds = slice_data[mask_matrix] - filter!(x -> x > 0, data_for_bounds) - else - data_for_bounds = vec(slice_data) - end - - # Handle cases where data_for_bounds might be empty after filtering - if isempty(data_for_bounds) - min_val, max_val = 0.0, 1.0 # Default range if no data in ROI - else - min_val, max_val = if use_triq - get_outlier_thres(data_for_bounds, triq_prob) +function generate_colorbar_image(slice_data::AbstractMatrix, color_levels::Int, output_path::String, + bounds::Tuple{Float64, Float64}; + use_triq::Bool=false, triq_prob::Float64=0.98, + mask_path::Union{String, Nothing}=nothing) + # Use the provided bounds instead of recalculating + min_val, max_val = bounds + + # If bounds are invalid, calculate fallback + if min_val == max_val == 0.0 + local data_for_bounds + if mask_path !== nothing + height, width = size(slice_data) + mask_matrix = load_and_prepare_mask(mask_path, (width, height)) + data_for_bounds = slice_data[mask_matrix] + filter!(x -> x > 0, data_for_bounds) else - extrema(data_for_bounds) + data_for_bounds = vec(slice_data) end + min_val, max_val = isempty(data_for_bounds) ? (0.0, 1.0) : extrema(data_for_bounds) end # 2. Replicate the tick calculation logic from plot_slices diff --git a/src/mzML.jl b/src/mzML.jl index aac75e9..1709780 100644 --- a/src/mzML.jl +++ b/src/mzML.jl @@ -36,10 +36,13 @@ determine the data type, compression, and axis type. function get_spectrum_asset_metadata(stream::IO) start_pos = position(stream) - bda_tag = find_tag(stream, r"(\d+)", footer) if index_offset_match === nothing - error("Could not find . File may not be an indexed mzML.") + throw(FileFormatError("Could not find . File may not be an indexed mzML.")) end return parse(Int64, index_offset_match.captures[1]) @@ -235,14 +238,14 @@ function load_mzml_lazy(file_path::String; cache_size::Int=100) println("DEBUG: Searching for ''.") if find_tag(stream, r"