diff --git a/julia_imzML_visual.jl b/julia_imzML_visual.jl
index 439cd14..aff0afa 100644
--- a/julia_imzML_visual.jl
+++ b/julia_imzML_visual.jl
@@ -626,4 +626,151 @@ function warmup_init()
println("Pre-compilation finished.")
end
+end
+
+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}()
+ end
+ end
+ return Dict{String, Any}()
+end
+
+function extract_metadata(msi_data::MSIData, source_path::String)
+ df = msi_data.spectrum_stats_df
+ if df === nothing
+ # This can happen if precompute_analytics hasn't been run
+ # We can still return basic info
+ return Dict(
+ "summary" => [
+ Dict("parameter" => "File Name", "value" => basename(source_path)),
+ Dict("parameter" => "Number of Spectra", "value" => length(msi_data.spectra_metadata)),
+ Dict("parameter" => "Image Dimensions", "value" => "$(msi_data.image_dims[1]) x $(msi_data.image_dims[2])"),
+ ],
+ "global_min_mz" => nothing,
+ "global_max_mz" => nothing
+ )
+ end
+
+ summary_stats = [
+ Dict("parameter" => "File Name", "value" => basename(source_path)),
+ Dict("parameter" => "Number of Spectra", "value" => length(msi_data.spectra_metadata)),
+ Dict("parameter" => "Image Dimensions", "value" => "$(msi_data.image_dims[1]) x $(msi_data.image_dims[2])"),
+ Dict("parameter" => "Global Min m/z", "value" => @sprintf("%.4f", msi_data.global_min_mz)),
+ Dict("parameter" => "Global Max m/z", "value" => @sprintf("%.4f", msi_data.global_max_mz)),
+ Dict("parameter" => "Mean TIC", "value" => @sprintf("%.2e", mean(df.TIC))),
+ Dict("parameter" => "Mean BPI", "value" => @sprintf("%.2e", mean(df.BPI))),
+ Dict("parameter" => "Mean # Points", "value" => @sprintf("%.1f", mean(df.NumPoints))),
+ ]
+
+ if hasproperty(df, :Mode)
+ centroid_count = count(==(MSI_src.CENTROID), df.Mode)
+ profile_count = count(==(MSI_src.PROFILE), df.Mode)
+ unknown_count = count(==(MSI_src.UNKNOWN), df.Mode)
+
+ push!(summary_stats, Dict("parameter" => "Centroid Spectra", "value" => string(centroid_count)))
+ push!(summary_stats, Dict("parameter" => "Profile Spectra", "value" => string(profile_count)))
+ if unknown_count > 0
+ push!(summary_stats, Dict("parameter" => "Unknown Mode Spectra", "value" => string(unknown_count)))
+ end
+ end
+
+ return Dict(
+ "summary" => summary_stats,
+ "global_min_mz" => msi_data.global_min_mz,
+ "global_max_mz" => msi_data.global_max_mz
+ )
+end
+
+function update_registry(registry_path, dataset_name, source_path, metadata=nothing, is_imzML=false)
+ registry = load_registry(registry_path)
+ entry = Dict{String, Any}( # Explicitly type the dictionary to allow mixed value types
+ "source_path" => source_path,
+ "processed_date" => string(now()),
+ "is_imzML" => is_imzML
+ )
+ if metadata !== nothing
+ entry["metadata"] = metadata
+ end
+ 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
+end
+
+function process_file_safely(file_path, masses, params, progress_message_ref, overall_progress_ref)
+ local_msi_data = nothing
+ dataset_name = replace(basename(file_path), r"\.imzML$"i => "")
+ output_dir = joinpath("public", dataset_name)
+ println("Processing: $dataset_name -> $output_dir")
+
+ try
+ # --- Load Data ---
+ progress_message_ref = "Loading: $(basename(file_path))"
+ local_msi_data = OpenMSIData(file_path)
+ if !(local_msi_data.source isa ImzMLSource)
+ @warn "Skipping non-imzML file: $(basename(file_path))"
+ return (false, "Skipped: Not an imzML file")
+ end
+
+ # --- Generate Slices (this will call precompute_analytics if needed) ---
+ progress_message_ref = "Generating $(length(masses)) slices for $(dataset_name)..."
+ slice_dict = get_multiple_mz_slices(local_msi_data, masses, params.tolerance)
+
+ # --- Extract metadata *after* it has been computed ---
+ metadata = extract_metadata(local_msi_data, file_path)
+
+ # --- Save Slices ---
+ mkpath(output_dir) # Ensure output directory exists
+
+ for (mass_idx, mass) in enumerate(masses)
+ progress_message_ref = "File $(params.fileIdx)/$(params.nFiles): Saving slice for m/z=$mass"
+
+ slice = slice_dict[mass]
+ text_nmass = replace(string(mass), "." => "_")
+ bitmap_filename = params.triqE ? "TrIQ_$(text_nmass).bmp" : "MSI_$(text_nmass).bmp"
+ colorbar_filename = params.triqE ? "colorbar_TrIQ_$(text_nmass).png" : "colorbar_MSI_$(text_nmass).png"
+
+ if all(iszero, slice)
+ sliceQuant = zeros(UInt8, size(slice))
+ @warn "No intensity data for m/z = $mass in $(dataset_name)"
+ else
+ sliceQuant = params.triqE ? TrIQ(slice, params.colorL, params.triqP) : quantize_intensity(slice, params.colorL)
+ if params.medianF
+ sliceQuant = round.(UInt8, median_filter(sliceQuant))
+ 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)
+ end
+ end
+
+ is_imzML = local_msi_data.source isa ImzMLSource
+ update_registry(params.registry, dataset_name, file_path, metadata, is_imzML)
+ return (true, "")
+
+ catch e
+ @error "File processing failed" file=file_path exception=(e, catch_backtrace())
+ return (false, "File: $(basename(file_path)) - $(sprint(showerror, e))")
+ finally
+ if local_msi_data !== nothing
+ # Cleanup
+ end
+ local_msi_data = nothing
+ GC.gc(true)
+ if Sys.islinux()
+ ccall(:malloc_trim, Int32, (Int32,), 0)
+ end
+ end
end
\ No newline at end of file
diff --git a/mask.jl b/mask.jl
new file mode 100644
index 0000000..349a5f2
--- /dev/null
+++ b/mask.jl
@@ -0,0 +1,897 @@
+module MaskApp
+
+# --- Dependencies ---
+using GenieFramework
+using Images, ImageBinarization, ImageMorphology, ImageComponentAnalysis
+using NativeFileDialog, FileIO, ImageCore, Printf, Dates, JSON
+using Pkg, Libz, PlotlyBase, CairoMakie, Colors, Base64
+using Statistics, NaturalSort, LinearAlgebra, StipplePlotly
+using Base.Filesystem: mv
+
+using MSI_src
+using .MSI_src: MSIData
+
+# Plot Handling
+include("./julia_imzML_visual.jl")
+
+# Image Processing Pipeline
+include("src/ImageProcessing.jl")
+using .ImageProcessing
+using ImageBinarization
+
+function load_and_binarize_mask(path)
+ img = load(path)
+ gray_img = ensure_grayscale(img)
+ # Binarize to ensure only pure black and white values
+ return binarize(gray_img, Otsu())
+end
+
+function ensure_grayscale(img)
+ if eltype(img) <: AbstractRGB; return Gray.(img); end
+ if eltype(img) <: AbstractRGBA; return Gray.(RGB.(img)); end
+ if eltype(img) <: Color; return Gray.(img); end
+ return img
+end
+
+function alter_image(img_path, otsu_scale, noise_size_percent, hole_size_percent, smoothing_level)
+ try
+ # Strip query parameters from img_path
+ clean_img_path = replace(img_path, r"\?.*" => "")
+ full_image_path = joinpath("public", lstrip(clean_img_path, '/'))
+ if !isfile(full_image_path)
+ return (success=false, message="Image file not found: $(full_image_path)", path="")
+ end
+
+ original_img = load(full_image_path)
+ gray_img = Float32.(ensure_grayscale(original_img))
+
+ binary, noise_removed, holes_filled, smoothed =
+ ImageProcessing.process_image_pipeline(gray_img;
+ otsu_scale=otsu_scale, noise_size_percent=noise_size_percent,
+ hole_size_percent=hole_size_percent, smoothing=smoothing_level)
+
+ output_dir = joinpath("public", "css", "masks")
+ mkpath(output_dir)
+
+ path_smooth = "/css/masks/smoothed.png"
+ save(joinpath(output_dir, basename(path_smooth)), smoothed)
+
+ return (success=true, message="Image processed successfully!", path=path_smooth)
+
+ catch e
+ @error "Mask editor processing failed" exception=(e, catch_backtrace())
+ return (success=false, message="Error processing image: $(sprint(showerror, e))", path="")
+ end
+end
+
+# New helper function to display a single slice
+function display_slice(slice_path)
+ if isempty(slice_path)
+ return [PlotlyBase.heatmap(x=Vector{Float64}(), y=Vector{Float64}())], PlotlyBase.Layout(margin=attr(l=0,r=0,t=0,b=0,pad=0)), false
+ end
+
+ full_path = joinpath("public", lstrip(slice_path, '/'))
+ if !isfile(full_path)
+ return [PlotlyBase.heatmap(x=Vector{Float64}(), y=Vector{Float64}())], PlotlyBase.Layout(margin=attr(l=0,r=0,t=0,b=0,pad=0)), false
+ end
+
+ plotdata, plotlayout, _, _ = loadImgPlot(slice_path, "", 0.0) # No mask overlay
+ return plotdata, plotlayout, true
+end
+
+function get_timestamped_path(base_path)
+ clean_path = replace(base_path, r"\?.*" => "")
+ return clean_path * "?t=" * string(time_ns())
+end
+
+function refresh_editor_preview(imgInt, smoothed_mask_path, imgTrans)
+ if !isempty(smoothed_mask_path)
+ timestamped_path = get_timestamped_path(smoothed_mask_path)
+ plotdata, plotlayout, _, _ = loadImgPlot(imgInt, timestamped_path, imgTrans)
+ return plotdata, plotlayout
+ else
+ return [PlotlyBase.heatmap(x=Vector{Float64}(), y=Vector{Float64}())], PlotlyBase.Layout(margin=attr(l=0,r=0,t=0,b=0,pad=0))
+ end
+end
+
+function flood_fill!(img::AbstractMatrix, x::Int, y::Int, fill_color)
+ h, w = size(img)
+ if !(1 <= x <= h && 1 <= y <= w); return; end
+
+ target_color = img[x, y]
+ if target_color == fill_color; return; end
+
+ q = [(x, y)]
+ img[x, y] = fill_color
+
+ while !isempty(q)
+ cx, cy = popfirst!(q)
+
+ for (dx, dy) in [(0, 1), (0, -1), (1, 0), (-1, 0)]
+ nx, ny = cx + dx, cy + dy
+ if 1 <= nx <= h && 1 <= ny <= w && img[nx, ny] == target_color
+ img[nx, ny] = fill_color
+ push!(q, (nx, ny))
+ end
+ end
+ end
+end
+
+function update_main_plot(imgInt::String, smoothed_mask_path::String, imgTrans::Float64)
+ try
+ if !isempty(imgInt) && !isempty(smoothed_mask_path)
+ # Both slice and mask available - show overlay
+ plotdata, plotlayout, _, _ = loadImgPlot(imgInt, smoothed_mask_path, imgTrans)
+ return plotdata, plotlayout, true
+ elseif !isempty(imgInt)
+ # Only slice available - show slice alone
+ plotdata, plotlayout, show = display_slice(imgInt)
+ return plotdata, plotlayout, show
+ elseif !isempty(smoothed_mask_path)
+ # Only mask available - show mask alone
+ plotdata, plotlayout, _, _ = loadImgPlot(smoothed_mask_path, "", 0.0)
+ return plotdata, plotlayout, true
+ else
+ # Nothing to show
+ return [PlotlyBase.heatmap(x=Vector{Float64}(), y=Vector{Float64}())],
+ PlotlyBase.Layout(margin=attr(l=0,r=0,t=0,b=0,pad=0)),
+ false
+ end
+ catch e
+ @error "Failed to update main plot" exception=(e, catch_backtrace())
+ return [PlotlyBase.heatmap(x=Vector{Float64}(), y=Vector{Float64}())],
+ PlotlyBase.Layout(margin=attr(l=0,r=0,t=0,b=0,pad=0)),
+ false
+ end
+end
+
+@genietools
+
+@app begin
+ # --- State for Slice Selection (from app.jl) ---
+ @in selected_folder_main = ""
+ @out available_folders = String[]
+ @out image_available_folders = String[]
+ @private registry_init_done = false
+
+ @out imgInt = "" # Path to current slice
+ @out current_msi = ""
+ @out msgimg = "Please select a dataset."
+ @in imgTrans = 0.5 # Default transparency
+
+ # --- State for Mask Workflow ---
+ @in otsu_scale = 1.0
+ @in noise_size_percent = 0.01
+ @in hole_size_percent = 0.005
+ @in smoothing_level = 3
+ @in is_editing_mask = false
+ @in is_browsing_slices = true
+
+ @in mask_input_path = "" # Reactive trigger for processing
+ @out smoothed_mask_path = "" # Path to the final smoothed.png
+
+ @out plotdata_verify = [PlotlyBase.heatmap(x=Vector{Float64}(), y=Vector{Float64}())]
+ @out plotlayout_verify = PlotlyBase.Layout(margin=attr(l=0,r=0,t=0,b=0,pad=0))
+ @out show_verification_plot = false
+
+ @out show_editor = false # To launch the custom editor
+
+ # --- Messages and Warnings ---
+ @out mask_editor_message = ""
+ @out mask_editor_warning = false
+ @out progress = false
+
+ # --- Buttons ---
+ @in btn_use_slice_as_mask = false
+ @in btn_upload_mask = false
+ @in btn_edit_manually = false
+ @in btn_img_plus = false
+ @in btn_img_minus = false
+ @in btn_change_slice = false # New button
+ @in btn_flip_mask = false
+ @in btn_save_final_mask = false
+ @in btn_move_mask = false
+
+ # --- Manual Editor State ---
+ @in brush_size = 10
+ @in brush_color = "#ff0000" # Red brush
+ @in current_tool = "brush" # "brush", "eraser", "bucket", "drag"
+ @in editor_scale = 1.0
+ @in bucket_fill_trigger = Dict()
+
+ # Canvas dimensions matching the smoothed mask
+ @out canvas_width = 512
+ @out canvas_height = 512
+
+ # Mouse coordinates for drawing
+ @in mouse_x = 0
+ @in mouse_y = 0
+ @in is_drawing = false
+
+ # Editor operations
+ @in rotate_degrees = 0
+ @in flip_direction = "horizontal"
+ @in move_direction = "right"
+ @in move_pixels = 10
+ @in move_mask_payload = Dict()
+
+ # Final save
+ @in final_mask_name = ""
+ @in updated_mask_data = "" # New property for canvas data
+ @in canvas_refresh_trigger = 0 # New property to trigger canvas refresh
+
+ # --- Handlers ---
+
+ @onchange updated_mask_data begin
+ if !isempty(updated_mask_data)
+ try
+ # Data URL format: data:image/png;base64,iVBORw0KGgo...
+ # Extract base64 part
+ base64_data = split(updated_mask_data, ",")[2]
+ decoded_img_bytes = base64decode(base64_data)
+
+ # Load image from bytes, convert to binary, and re-save
+ img_from_canvas = load(IOBuffer(decoded_img_bytes))
+ binary_mask = binarize(ensure_grayscale(img_from_canvas), Otsu())
+
+ mask_path = joinpath("public", lstrip(replace(smoothed_mask_path, r"\?.*" => ""), '/'))
+ save(mask_path, binary_mask)
+
+ smoothed_mask_path = get_timestamped_path(smoothed_mask_path)
+ @info "Mask updated from canvas data and binarized."
+ catch e
+ @error "Failed to update mask from canvas data" exception=(e, catch_backtrace())
+ end
+ end
+ end
+
+ @onchange selected_folder_main begin
+ if !isempty(selected_folder_main)
+ is_browsing_slices = true
+ is_editing_mask = false
+ folder_path = joinpath("public", selected_folder_main)
+ println("Selected folder: $selected_folder_main")
+ if !isdir(folder_path)
+ imgInt = ""
+ msgimg = "Folder not found."
+ show_verification_plot = false
+ return
+ end
+
+ msi_bmp = sort(filter(f -> startswith(f, "MSI_") && endswith(f, ".bmp"), readdir(folder_path)), lt=natural)
+
+ if !isempty(msi_bmp)
+ current_msi = first(msi_bmp)
+ imgInt = "/$(selected_folder_main)/$(current_msi)"
+ text_nmass = replace(current_msi, r"MSI_|.bmp" => "")
+ msgimg = "
m/z: $(replace(text_nmass, "_" => "."))"
+
+ # Use centralized plot update
+ plotdata_verify, plotlayout_verify, show_verification_plot = update_main_plot(imgInt, smoothed_mask_path, imgTrans)
+ else
+ imgInt = ""
+ msgimg = "No MSI images found in this dataset."
+ show_verification_plot = false
+ end
+ end
+ end
+
+ @onbutton btn_img_plus begin
+ if isempty(selected_folder_main) || is_editing_mask return end
+ folder_path = joinpath("public", selected_folder_main)
+ if !isdir(folder_path) return end
+ msi_bmp = sort(filter(f -> startswith(f, "MSI_") && endswith(f, ".bmp"), readdir(folder_path)), lt=natural)
+ new_msi = increment_image(current_msi, msi_bmp)
+ if new_msi !== nothing
+ current_msi = new_msi
+ imgInt = "/$(selected_folder_main)/$(current_msi)"
+ text_nmass = replace(current_msi, r"MSI_|.bmp" => "")
+ msgimg = "
m/z: $(replace(text_nmass, "_" => "."))"
+
+ # Use centralized plot update
+ plotdata_verify, plotlayout_verify, show_verification_plot = update_main_plot(imgInt, smoothed_mask_path, imgTrans)
+ end
+ end
+
+ @onbutton btn_img_minus begin
+ if isempty(selected_folder_main) || is_editing_mask return end
+ folder_path = joinpath("public", selected_folder_main)
+ if !isdir(folder_path) return end
+ msi_bmp = sort(filter(f -> startswith(f, "MSI_") && endswith(f, ".bmp"), readdir(folder_path)), lt=natural)
+ new_msi = decrement_image(current_msi, msi_bmp)
+ if new_msi !== nothing
+ current_msi = new_msi
+ imgInt = "/$(selected_folder_main)/$(current_msi)"
+ text_nmass = replace(current_msi, r"MSI_|.bmp" => "")
+ msgimg = "
m/z: $(replace(text_nmass, "_" => "."))"
+
+ # Use centralized plot update
+ plotdata_verify, plotlayout_verify, show_verification_plot = update_main_plot(imgInt, smoothed_mask_path, imgTrans)
+ end
+ end
+
+ @onbutton btn_use_slice_as_mask begin
+ if !isempty(imgInt)
+ is_browsing_slices = false
+ is_editing_mask = true
+ mask_input_path = imgInt
+ mask_editor_message = "Using current slice as mask input..."
+ mask_editor_warning = false
+
+ # Trigger processing which will update the plot via mask_input_path handler
+ else
+ mask_editor_message = "No slice selected. Please select a dataset and slice first."
+ mask_editor_warning = true
+ end
+ end
+
+ @onbutton btn_upload_mask begin
+ picked_path = pick_file(; filterlist="png,bmp,jpg,jpeg")
+ if !isempty(picked_path)
+ target_dir = joinpath("public", "css", "masks")
+ mkpath(target_dir)
+ uploaded_filename = "uploaded_mask_input.png"
+ destination_path = joinpath(target_dir, uploaded_filename)
+ cp(picked_path, destination_path; force=true)
+
+ is_browsing_slices = false
+ is_editing_mask = true
+ mask_input_path = get_timestamped_path("/css/masks/$uploaded_filename")
+ mask_editor_message = "Uploaded mask image for processing..."
+ mask_editor_warning = false
+
+ # Plot will be updated via mask_input_path handler after processing
+ else
+ mask_editor_message = "No file selected for upload."
+ mask_editor_warning = true
+ end
+ end
+
+ @onbutton btn_change_slice begin
+ is_browsing_slices = true
+ is_editing_mask = false
+
+ # Use centralized plot update
+ plotdata_verify, plotlayout_verify, show_verification_plot = update_main_plot(imgInt, "", imgTrans)
+ smoothed_mask_path = "" # Clear the mask path
+ mask_editor_message = "Slice browsing re-enabled."
+ end
+
+ @onchange mask_input_path begin
+ if isempty(mask_input_path) return end
+
+ result = alter_image(mask_input_path, otsu_scale, noise_size_percent, hole_size_percent, smoothing_level)
+
+ if result.success
+ local_smoothed_path = joinpath("public", lstrip(result.path, '/'))
+ slice_path_cleaned = replace(imgInt, r"\?.*" => "")
+ slice_full_path = joinpath("public", lstrip(slice_path_cleaned, '/'))
+
+ if isfile(slice_full_path) && isfile(local_smoothed_path)
+ slice_img = load(slice_full_path)
+ mask_img = load(local_smoothed_path)
+
+ if size(slice_img) != size(mask_img)
+ @info "Resizing mask to match slice dimensions."
+ resized_mask = imresize(mask_img, size(slice_img))
+ save(local_smoothed_path, resized_mask)
+ end
+
+ smoothed_mask_path = get_timestamped_path(result.path)
+
+ # Use centralized plot update
+ plotdata_verify, plotlayout_verify, show_verification_plot = update_main_plot(imgInt, smoothed_mask_path, imgTrans)
+ mask_editor_message = result.message
+ mask_editor_warning = false
+ elseif isfile(local_smoothed_path)
+ smoothed_mask_path = get_timestamped_path(result.path)
+ # Use centralized plot update
+ plotdata_verify, plotlayout_verify, show_verification_plot = update_main_plot("", smoothed_mask_path, imgTrans)
+ mask_editor_message = "Displaying generated mask. Select a dataset to overlay a slice."
+ mask_editor_warning = false
+ else
+ mask_editor_message = "Slice or mask image not found for verification plot."
+ mask_editor_warning = true
+ end
+ else
+ mask_editor_message = result.message
+ mask_editor_warning = true
+ end
+ end
+
+ @onchange otsu_scale, noise_size_percent, hole_size_percent, smoothing_level begin
+ if is_editing_mask
+ if isempty(mask_input_path) return end
+
+ result = alter_image(mask_input_path, otsu_scale, noise_size_percent, hole_size_percent, smoothing_level)
+
+ if result.success
+ local_smoothed_path = joinpath("public", lstrip(result.path, '/'))
+ slice_path_cleaned = replace(imgInt, r"\?.*" => "")
+ slice_full_path = joinpath("public", lstrip(slice_path_cleaned, '/'))
+
+ if isfile(slice_full_path) && isfile(local_smoothed_path)
+ slice_img = load(slice_full_path)
+ mask_img = load(local_smoothed_path)
+
+ if size(slice_img) != size(mask_img)
+ @info "Resizing mask to match slice dimensions."
+ resized_mask = imresize(mask_img, size(slice_img))
+ save(local_smoothed_path, resized_mask)
+ end
+
+ smoothed_mask_path = get_timestamped_path(result.path)
+
+ # Use centralized plot update
+ plotdata_verify, plotlayout_verify, show_verification_plot = update_main_plot(imgInt, smoothed_mask_path, imgTrans)
+ mask_editor_message = "Mask updated."
+ mask_editor_warning = false
+ elseif isfile(local_smoothed_path)
+ smoothed_mask_path = get_timestamped_path(result.path)
+ # Use centralized plot update
+ plotdata_verify, plotlayout_verify, show_verification_plot = update_main_plot("", smoothed_mask_path, imgTrans)
+ mask_editor_message = "Displaying generated mask. Select a dataset to overlay a slice."
+ mask_editor_warning = false
+ else
+ mask_editor_message = "Slice or mask image not found for verification plot."
+ mask_editor_warning = true
+ end
+ else
+ mask_editor_message = result.message
+ mask_editor_warning = true
+ end
+ end
+ end
+
+ @onchange imgTrans begin
+ if is_editing_mask
+ # Use centralized plot update
+ plotdata_verify, plotlayout_verify, show_verification_plot = update_main_plot(imgInt, smoothed_mask_path, imgTrans)
+ end
+ end
+
+ @onbutton btn_edit_manually begin
+ if !isempty(smoothed_mask_path)
+ mask_path = joinpath("public", lstrip(replace(smoothed_mask_path, r"\?.*" => ""), '/'))
+ if isfile(mask_path)
+ mask_img = load(mask_path)
+ h, w = size(mask_img)
+ canvas_width = w
+ canvas_height = h
+ show_editor = true
+ end
+ else
+ mask_editor_message = "No mask available for editing. Process a mask first."
+ mask_editor_warning = true
+ end
+ end
+
+ @onchange show_editor begin
+ if !show_editor && is_editing_mask
+ # When the editor dialog is closed, refresh the main plot using centralized function
+ plotdata_verify, plotlayout_verify, show_verification_plot = update_main_plot(imgInt, smoothed_mask_path, imgTrans)
+ is_drawing = false # Reset drawing state
+ end
+ end
+
+ @onchange mouse_x, mouse_y, is_drawing begin
+ if show_editor && is_drawing && mouse_x > 0 && mouse_y > 0
+ try
+ mask_path = joinpath("public", lstrip(replace(smoothed_mask_path, r"\?.*" => ""), '/'))
+ if !isfile(mask_path) return end
+
+ img = load(mask_path)
+
+ actual_x = round(Int, mouse_x)
+ actual_y = round(Int, mouse_y)
+
+ if current_tool == "brush" || current_tool == "eraser"
+ if brush_size <= 0 return end
+ radius = brush_size / 2 # Treat as diameter
+
+ for i in max(1, floor(Int, actual_y-radius)):min(size(img, 1), ceil(Int, actual_y+radius))
+ for j in max(1, floor(Int, actual_x-radius)):min(size(img, 2), ceil(Int, actual_x+radius))
+ if (i - actual_y)^2 + (j - actual_x)^2 <= radius^2
+ color = (current_tool == "brush") ? RGB(1, 1, 1) : RGB(0, 0, 0)
+ img[i, j] = color
+ end
+ end
+ end
+ end
+
+ save(mask_path, img)
+ plotdata_verify, plotlayout_verify, show_verification_plot = update_main_plot(imgInt, smoothed_mask_path, imgTrans)
+ smoothed_mask_path = get_timestamped_path(smoothed_mask_path)
+
+ catch e
+ @error "Editor action failed" exception=(e, catch_backtrace())
+ end
+ end
+ end
+
+ @onchange bucket_fill_trigger begin
+ if show_editor && !isempty(bucket_fill_trigger)
+ try
+ mask_path = joinpath("public", lstrip(replace(smoothed_mask_path, r"\?.*" => ""), '/'))
+ if !isfile(mask_path) return end
+
+ img = load_and_binarize_mask(mask_path)
+
+ actual_x = round(Int, bucket_fill_trigger["x"])
+ actual_y = round(Int, bucket_fill_trigger["y"])
+
+ if 1 <= actual_y <= size(img, 1) && 1 <= actual_x <= size(img, 2)
+ fill_color = Gray(1) # White
+ flood_fill!(img, actual_y, actual_x, fill_color)
+
+ save(mask_path, img)
+ plotdata_verify, plotlayout_verify, show_verification_plot = update_main_plot(imgInt, smoothed_mask_path, imgTrans)
+ smoothed_mask_path = get_timestamped_path(smoothed_mask_path)
+ canvas_refresh_trigger = canvas_refresh_trigger[] + 1
+ end
+ catch e
+ @error "Bucket fill failed" exception=(e, catch_backtrace())
+ end
+ end
+ end
+
+ @onchange rotate_degrees begin
+ if show_editor && rotate_degrees != 0
+ try
+ mask_path = joinpath("public", lstrip(replace(smoothed_mask_path, r"\?.*" => ""), '/'))
+ if !isfile(mask_path) return end
+
+ img = load_and_binarize_mask(mask_path)
+ original_size = size(img)
+
+ # Normalize angle to 0, 90, 180, 270
+ angle = mod(round(Int, rotate_degrees), 360)
+ if angle == 0; rotate_degrees = 0; return; end
+
+ rotated_img = if angle == 90
+ rot_r90(img)
+ elseif angle == 180
+ rot180(img)
+ elseif angle == 270
+ rot_l90(img)
+ else
+ img # No change for other angles
+ end
+
+ if rotated_img !== img
+ # Create a new image with the original dimensions, filled with black
+ new_img = similar(img, original_size)
+ fill!(new_img, eltype(img)(0))
+
+ # Calculate padding/offset to center the rotated image
+ rotated_size = size(rotated_img)
+ offset_h = (original_size[1] - rotated_size[1]) ÷ 2
+ offset_w = (original_size[2] - rotated_size[2]) ÷ 2
+
+ # Define the region in the new image where the rotated image will be placed
+ dest_region_h = (1:rotated_size[1]) .+ offset_h
+ dest_region_w = (1:rotated_size[2]) .+ offset_w
+
+ # Ensure the destination region is within the bounds of the new image
+ clamped_dest_h = max(1, dest_region_h.start):min(original_size[1], dest_region_h.stop)
+ clamped_dest_w = max(1, dest_region_w.start):min(original_size[2], dest_region_w.stop)
+
+ # Define the source region from the rotated image
+ src_h = (1:length(clamped_dest_h))
+ src_w = (1:length(clamped_dest_w))
+
+ if !isempty(clamped_dest_h) && !isempty(clamped_dest_w)
+ view(new_img, clamped_dest_h, clamped_dest_w) .= view(rotated_img, src_h, src_w)
+ end
+
+ save(mask_path, new_img)
+ plotdata_verify, plotlayout_verify, show_verification_plot = update_main_plot(imgInt, smoothed_mask_path, imgTrans)
+ end
+
+ # Reset the slider and update the path
+ smoothed_mask_path = get_timestamped_path(smoothed_mask_path)
+ rotate_degrees = 0
+ canvas_refresh_trigger = canvas_refresh_trigger[] + 1
+
+ catch e
+ @error "Rotation failed" exception=(e, catch_backtrace())
+ end
+ end
+ end
+
+ @onbutton btn_flip_mask begin
+ if show_editor
+ try
+ mask_path = joinpath("public", lstrip(replace(smoothed_mask_path, r"\?.*" => ""), '/'))
+ if !isfile(mask_path) return end
+
+ img = load_and_binarize_mask(mask_path)
+ flipped = (flip_direction == "horizontal") ? reverse(img, dims=2) : reverse(img, dims=1)
+ save(mask_path, flipped)
+ plotdata_verify, plotlayout_verify, show_verification_plot = update_main_plot(imgInt, smoothed_mask_path, imgTrans)
+
+ smoothed_mask_path = get_timestamped_path(smoothed_mask_path)
+ canvas_refresh_trigger = canvas_refresh_trigger[] + 1
+
+ catch e
+ @error "Flip failed" exception=(e, catch_backtrace())
+ end
+ end
+ end
+
+ @onchange move_mask_payload begin
+ if show_editor && move_pixels > 0 && !isempty(move_mask_payload)
+ try
+ mask_path = joinpath("public", lstrip(replace(smoothed_mask_path, r"\?.*" => ""), '/'))
+ if !isfile(mask_path) return end
+
+ img = load_and_binarize_mask(mask_path)
+ h, w = size(img)
+ new_img = similar(img)
+ fill!(new_img, eltype(img)(0))
+
+ dx, dy = 0, 0
+ direction = move_mask_payload["direction"]
+ if direction == "right"; dx = move_pixels; end
+ if direction == "left"; dx = -move_pixels; end
+ if direction == "down"; dy = move_pixels; end
+ if direction == "up"; dy = -move_pixels; end
+
+ src_y_range = max(1, 1-dy):min(h, h-dy)
+ src_x_range = max(1, 1-dx):min(w, w-dx)
+ dest_y_range = max(1, 1+dy):min(h, h+dy)
+ dest_x_range = max(1, 1+dx):min(w, w+dx)
+
+ if !isempty(src_y_range) && !isempty(src_x_range) && !isempty(dest_y_range) && !isempty(dest_x_range)
+ view(new_img, dest_y_range, dest_x_range) .= view(img, src_y_range, src_x_range)
+ end
+
+ save(mask_path, new_img)
+ plotdata_verify, plotlayout_verify, show_verification_plot = update_main_plot(imgInt, smoothed_mask_path, imgTrans)
+ smoothed_mask_path = get_timestamped_path(smoothed_mask_path)
+ canvas_refresh_trigger = canvas_refresh_trigger[] + 1
+
+ catch e
+ @error "Move failed" exception=(e, catch_backtrace())
+ end
+ end
+ end
+
+ @onbutton btn_save_final_mask begin
+ try
+ parent_folder = split(selected_folder_main, '/')[end]
+ final_mask_name = "$(parent_folder).png"
+
+ source_path = joinpath("public", lstrip(replace(smoothed_mask_path, r"\?.*" => ""), '/'))
+ target_dir = joinpath("public", "css", "masks")
+ mkpath(target_dir)
+ final_path = joinpath(target_dir, final_mask_name)
+
+ if isfile(source_path)
+ # Copy the mask to final location
+ cp(source_path, final_path; force=true)
+
+ # Update registry directly in the handler
+ reg_path = abspath(joinpath(@__DIR__, "public", "registry.json"))
+ registry = isfile(reg_path) ? JSON.parsefile(reg_path) : Dict{String, Any}()
+
+ # Update the registry entry
+ if haskey(registry, selected_folder_main)
+ registry[selected_folder_main]["mask_path"] = "/css/masks/$(final_mask_name)"
+ registry[selected_folder_main]["has_mask"] = true
+ else
+ # Create a new entry if folder doesn't exist in registry
+ registry[selected_folder_main] = Dict(
+ "mask_path" => "/css/masks/$(final_mask_name)",
+ "has_mask" => true,
+ "is_imzML" => true,
+ "processed_date" => string(Dates.now())
+ )
+ end
+
+ # Save the updated registry
+ open(reg_path, "w") do f
+ JSON.print(f, registry, 4)
+ end
+
+ @info "Registry updated with mask: $(final_mask_name)"
+
+ # Update UI state
+ mask_editor_message = "Final mask saved as $(final_mask_name)"
+ mask_editor_warning = false
+ show_editor = false
+ else
+ mask_editor_message = "Source mask file not found: $(source_path)"
+ mask_editor_warning = true
+ end
+
+ catch e
+ @error "Save final mask failed" exception=(e, catch_backtrace())
+ mask_editor_message = "Error saving final mask: $(sprint(showerror, e))"
+ mask_editor_warning = true
+ end
+ 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}()
+
+ 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
+ end
+ end
+ end
+ end
+
+ @methods """
+ // Close the default methods object and define our own component structure
+ },
+
+ created() {
+ // Initialize non-reactive data properties here
+ this.canvas = null;
+ this.ctx = null;
+ this.isDrawingOnCanvas = false;
+ this.lastX = 0;
+ this.lastY = 0;
+ this.imgObj = null;
+ },
+
+ watch: {
+ show_editor(newValue) {
+ if (newValue) {
+ // Wait for the dialog to render before initializing canvas
+ setTimeout(() => { this.initCanvas() }, 100);
+ }
+ },
+ canvas_refresh_trigger(newValue) {
+ if (newValue > 0) {
+ this.initCanvas();
+ }
+ }
+ },
+
+ methods: {
+ // Re-opened methods object for all our functions
+ initCanvas() {
+ this.canvas = document.getElementById('maskCanvas');
+ if (!this.canvas) {
+ console.error("Canvas element not found!");
+ return;
+ }
+ this.ctx = this.canvas.getContext('2d');
+
+ this.ctx.clearRect(0, 0, this.canvas.width, this.canvas.height);
+
+ this.imgObj = new Image();
+ this.imgObj.onload = () => {
+ this.ctx.drawImage(this.imgObj, 0, 0, this.canvas.width, this.canvas.height);
+ };
+ this.imgObj.src = this.smoothed_mask_path.split('?')[0] + '?t=' + Date.now();
+ },
+
+ startDrawing(event) {
+ if (this.current_tool === 'drag' || !this.ctx) return;
+
+ this.isDrawingOnCanvas = true;
+ const { x, y } = this.getCanvasMousePosition(event);
+ this.lastX = x;
+ this.lastY = y;
+
+ if (this.current_tool === 'bucket') {
+ this.bucket_fill_trigger = { x: x, y: y, t: Date.now() };
+ this.isDrawingOnCanvas = false;
+ return;
+ }
+
+ this.ctx.beginPath();
+ this.ctx.moveTo(this.lastX, this.lastY);
+ },
+
+ draw(event) {
+ if (!this.isDrawingOnCanvas || this.current_tool === 'drag' || this.current_tool === 'bucket' || !this.ctx) return;
+
+ const { x, y } = this.getCanvasMousePosition(event);
+
+ this.ctx.lineWidth = this.brush_size;
+ this.ctx.lineCap = 'round';
+ this.ctx.lineJoin = 'round';
+
+ if (this.current_tool === 'brush') {
+ this.ctx.strokeStyle = 'white';
+ this.ctx.globalCompositeOperation = 'source-over';
+ } else if (this.current_tool === 'eraser') {
+ this.ctx.strokeStyle = 'black';
+ this.ctx.globalCompositeOperation = 'source-over';
+ }
+
+ this.ctx.lineTo(x, y);
+ this.ctx.stroke();
+ },
+
+ stopDrawing() {
+ if (this.isDrawingOnCanvas) {
+ this.isDrawingOnCanvas = false;
+ this.ctx.closePath();
+ this.syncMaskToServer();
+ }
+ },
+
+ getCanvasMousePosition(event) {
+ const rect = this.canvas.getBoundingClientRect();
+ const scaleX = this.canvas.width / rect.width;
+ const scaleY = this.canvas.height / rect.height;
+ const x = (event.clientX - rect.left) * scaleX;
+ const y = (event.clientY - rect.top) * scaleY;
+ return { x, y };
+ },
+
+ syncMaskToServer() {
+ if (!this.canvas) return;
+ const dataURL = this.canvas.toDataURL('image/png');
+ this.updated_mask_data = dataURL;
+ },
+
+ // Kept for compatibility with bucket tool logic which uses original image coordinates
+ updateMousePosition(event) {
+ const rect = event.target.getBoundingClientRect();
+ const scale = this.editor_scale || 1;
+ this.mouse_x = (event.clientX - rect.left) / scale;
+ this.mouse_y = (event.clientY - rect.top) / scale;
+ }
+ // The closing brace for methods is intentionally omitted, as Genie/Stipple will add it.
+ """
+
+end
+
+@page("/mask", "mask.jl.html")
+
+end # module
diff --git a/mask.jl.html b/mask.jl.html
new file mode 100644
index 0000000..5d5492d
--- /dev/null
+++ b/mask.jl.html
@@ -0,0 +1,145 @@
+
+
+
+
+
+
+
Mask Generation Workflow
+
+
+
Step 1: Select Slice
+
+
+
+
+
+
+
Step 2: Create/Upload Mask
+
+
+
+
+
+
+
+
Step 3: Adjust & Verify
+
+
Adjust automated processing for the initial mask:
+
+
+
+
+
+
+
Adjust slice overlay transparency:
+
+
+
+
{{ mask_editor_message }}
+
+
+
+
+
+
+
+
+
+
+
Mask Preview
+
+
+
+
+
+
+
+
Please select a dataset and provide a mask input to begin.
+
Use the controls on the left to load a slice or upload a mask image.
+
+
+
+
+
+
+
+
+
+
+
+ Manual Mask Editor
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
diff --git a/public/css/genieapp.css b/public/css/genieapp.css
index f50c1f7..98ddaf2 100644
--- a/public/css/genieapp.css
+++ b/public/css/genieapp.css
@@ -90,6 +90,12 @@
overflow-y: auto; /* Add scroll if content overflows */
}
+.pixelated-plot svg image {
+ image-rendering: pixelated;
+ image-rendering: -moz-crisp-edges; /* Firefox */
+ image-rendering: crisp-edges;
+}
+
*{
font-family: 'Roboto', 'Lato', sans-serif;
}
\ No newline at end of file
diff --git a/routes.jl b/routes.jl
new file mode 100644
index 0000000..0d7c32c
--- /dev/null
+++ b/routes.jl
@@ -0,0 +1,3 @@
+# Include all your app modules
+include("app.jl")
+include("mask.jl")
diff --git a/src/ImageProcessing.jl b/src/ImageProcessing.jl
new file mode 100644
index 0000000..2f0f8a9
--- /dev/null
+++ b/src/ImageProcessing.jl
@@ -0,0 +1,47 @@
+module ImageProcessing
+
+using Images
+using ImageBinarization
+using ImageMorphology
+using ImageComponentAnalysis
+
+export process_image_pipeline
+
+# ===================================================================
+# CORE PROCESSING PIPELINE
+# ===================================================================
+
+function process_image_pipeline(gray_img;
+ otsu_scale=1.0,
+ noise_size_percent=0.1,
+ hole_size_percent=0.05,
+ smoothing=2)
+
+ # --- Step 1: Otsu Binarization ---
+ otsu_threshold = find_threshold(gray_img, Otsu())
+ adjusted_threshold = otsu_threshold * otsu_scale
+ binary_img = gray_img .>= adjusted_threshold
+
+ # --- Smart Parameter Scaling ---
+ image_area = length(gray_img)
+ noise_size_pixels = round(Int, image_area * noise_size_percent)
+ hole_size_pixels = round(Int, image_area * hole_size_percent)
+
+ # --- Step 2: Remove Small White Regions (Noise) ---
+ # area_opening is the correct morphological operation for this.
+ noise_removed_img = area_opening(binary_img, min_area=noise_size_pixels)
+
+ # --- Step 3: Fill Small Black Holes ---
+ # area_closing is the dual of area_opening and fills holes.
+ holes_filled_img = area_closing(noise_removed_img, min_area=hole_size_pixels)
+
+ # --- Step 4: Smooth Edges ---
+ # A morphological closing with a small disk smooths outlines.
+ smoothing_kernel = ones(Bool, (smoothing, smoothing))
+ smoothed_img = closing(holes_filled_img, smoothing_kernel)
+
+ # --- Return all intermediate steps for visualization ---
+ return binary_img, noise_removed_img, holes_filled_img, smoothed_img
+end
+
+end
diff --git a/src/MSI_src.jl b/src/MSI_src.jl
index 444dce5..19c0cf6 100644
--- a/src/MSI_src.jl
+++ b/src/MSI_src.jl
@@ -1,7 +1,40 @@
+# MSI_src.jl
+
module MSI_src
-# Export the public API
-export OpenMSIData, GetSpectrum, IterateSpectra, ImportMzmlFile, load_slices, plot_slices, plot_slice, get_total_spectrum, get_average_spectrum, LoadMzml, LoadSpectra, precompute_analytics
+# Export the public MSI API
+export OpenMSIData,
+ GetSpectrum,
+ IterateSpectra,
+ ImportMzmlFile,
+ load_slices,
+ plot_slices,
+ plot_slice,
+ get_total_spectrum,
+ get_average_spectrum,
+ LoadMzml,
+ precompute_analytics,
+ process_spectrum
+
+# Export the public Preprocessing API
+export FeatureMatrix,
+ run_preprocessing_pipeline,
+ qc_is_empty,
+ qc_is_regular,
+ transform_intensity,
+ smooth_spectrum,
+ snip_baseline,
+ tic_normalize,
+ pqn_normalize,
+ detect_peaks_profile,
+ align_peaks_lowess,
+ bin_peaks,
+ plot_stage_spectrum,
+ calculate_ppm_error,
+ calculate_resolution_fwhm,
+ analyze_mass_accuracy,
+ generate_qc_report,
+ get_common_calibration_standards
# Include all source files directly into the main module
include("MSIData.jl")
@@ -9,6 +42,7 @@ include("ParserHelpers.jl")
include("mzML.jl")
include("imzML.jl")
include("MzmlConverter.jl")
+include("Preprocessing.jl")
# --- Main Entry Point --- #
diff --git a/src/MzmlConverter.jl b/src/MzmlConverter.jl
index 4142c41..129cdf3 100644
--- a/src/MzmlConverter.jl
+++ b/src/MzmlConverter.jl
@@ -1,3 +1,4 @@
+# src/MzmlConverter.jl
"""
This file contains the workflow for converting .mzML files (with one spectrum per pixel)
into a proper .imzML/.ibd file pair, using a separate synchronization file.
diff --git a/src/Preprocessing.jl b/src/Preprocessing.jl
new file mode 100644
index 0000000..62e1e30
--- /dev/null
+++ b/src/Preprocessing.jl
@@ -0,0 +1,979 @@
+# src/Preprocessing.jl
+"""
+This module provides a comprehensive workflow for mass spectrometry imaging (MSI) data
+preprocessing, inspired by the functionality of the R package MALDIquant. It includes
+functions for quality control, intensity transformation, smoothing, baseline correction,
+normalization, peak picking, alignment, and feature matrix generation.
+"""
+
+# =============================================================================
+# Dependencies
+# =============================================================================
+
+using Statistics # For mean, median
+using StatsBase # For mad (Median Absolute Deviation)
+using SavitzkyGolay # For SavitzkyGolay filtering
+using Dates # For now()
+using CSV # For writing CSV files
+using DataFrames # For creating dataframes
+
+# =============================================================================
+# Data Structures
+# =============================================================================
+
+"""
+ FeatureMatrix
+
+Structure to hold the final feature matrix, including m/z bin boundaries and sample indices.
+
+# Fields
+- `matrix`: The numerical matrix where rows are samples and columns are features (m/z bins).
+- `mz_bins`: A vector of tuples `(low_mz, high_mz)` for each feature column.
+- `sample_ids`: A vector of indices corresponding to the original spectra.
+"""
+struct FeatureMatrix
+ matrix::Array{Float64,2} # samples × features
+ mz_bins::Vector{Tuple{Float64,Float64}} # [(low, high), ...]
+ sample_ids::Vector{Int}
+end
+
+# =============================================================================
+# 0) Quality Control (QC)
+# =============================================================================
+
+"""
+ qc_is_empty(mz, intensity) -> Bool
+
+Returns `true` if the spectrum is empty or contains no finite intensity data.
+"""
+qc_is_empty(mz::AbstractVector, intensity::AbstractVector)::Bool =
+ isempty(mz) || isempty(intensity) || all(!isfinite, intensity)
+
+"""
+ qc_is_regular(mz) -> Bool
+
+Checks that the m/z axis is monotonically non-decreasing, as expected in a profile spectrum.
+"""
+function qc_is_regular(mz::AbstractVector)
+ n = length(mz)
+ n < 2 && return true
+ @inbounds for i in 2:n
+ if mz[i] < mz[i-1]
+ return false
+ end
+ end
+ return true
+end
+
+# =============================================================================
+# 1) Intensity Transformation & Smoothing
+# =============================================================================
+
+"""
+ transform_intensity(intensity; method=:sqrt) -> Vector
+
+Applies a variance-stabilizing transformation to the intensity vector.
+Supported methods: `:sqrt` (default) and `:log1p`.
+"""
+function transform_intensity(intensity::AbstractVector{<:Real}; method::Symbol=:sqrt)
+ if method === :sqrt
+ return sqrt.(max.(zero(eltype(intensity)), intensity))
+ elseif method === :log1p
+ return log1p.(max.(zero(eltype(intensity)), intensity))
+ else
+ return collect(float.(intensity))
+ end
+end
+
+"""
+ smooth_spectrum(y; window=21, order=2) -> Vector
+
+Applies a Savitzky–Golay filter if `SavitzkyGolay.jl` is available.
+Otherwise, falls back to a simple (non-phase-correct) moving average.
+"""
+function smooth_spectrum(y::AbstractVector{<:Real}; window::Int=9, order::Int=2)
+ win = isodd(window) ? window : window + 1
+ res = SavitzkyGolay.savitzky_golay(collect(float.(y)), win, order)
+ return res.y
+end
+
+# =============================================================================
+# 2) Baseline Correction
+# =============================================================================
+
+"""
+ snip_baseline(y, iterations=100) -> Vector
+
+Estimates the baseline using a simple 1D SNIP (Statistics-sensitive Non-linear
+Iterative Peak-clipping) algorithm. `iterations` controls the aggressiveness.
+"""
+function snip_baseline(y::AbstractVector{<:Real}, iterations::Int=100)
+ n = length(y)
+ b = collect(float.(y)) # work copy
+ buf = similar(b)
+ for k in 1:iterations
+ copyto!(buf, b)
+ @inbounds for i in 2:n-1
+ buf[i] = min(b[i], 0.5 * (b[i-1] + b[i+1]))
+ end
+ # Handle endpoints
+ buf[1] = min(b[1], b[2])
+ buf[end] = min(b[end], b[end-1])
+ b, buf = buf, b # Swap buffers
+ end
+ return b
+end
+
+# =============================================================================
+# 3) Intensity Normalization
+# =============================================================================
+
+"""
+ tic_normalize(y) -> Vector
+
+Normalizes intensities to the Total Ion Current (TIC). If the sum is zero, returns a copy.
+"""
+function tic_normalize(y::AbstractVector{<:Real})
+ s = sum(y)
+ return s <= 0 ? collect(float.(y)) : collect(float.(y)) ./ s
+end
+
+"""
+ pqn_normalize(M) -> Matrix
+
+Performs Probabilistic Quotient Normalization on a matrix `M` where columns are spectra.
+"""
+function pqn_normalize(M::AbstractMatrix{<:Real})
+ M_float = collect(float.(M))
+ # Calculate reference spectrum (median across all spectra)
+ ref = mapslices(median, M_float; dims=2)[:,1]
+
+ # Calculate quotients for each spectrum relative to the reference
+ Q = similar(M_float)
+ @inbounds for j in axes(M_float, 2)
+ Q[:, j] = M_float[:, j] ./ (ref .+ eps(eltype(M_float)))
+ end
+
+ # Find the median quotient for each spectrum (scaling factor)
+ s = [median( @view Q[:, j]) for j in axes(Q, 2)]
+
+ # Normalize the original matrix
+ @inbounds for j in axes(M_float, 2)
+ M_float[:, j] ./= (s[j] + eps(eltype(M_float)))
+ end
+ return M_float
+end
+
+# =============================================================================
+# 4) Peak Detection (for Profile Data)
+# =============================================================================
+
+"""
+ detect_peaks_profile(mz, y; half_window=10, snr_threshold=2.0)
+
+Detects local maxima with a signal-to-noise threshold (using MAD for noise estimation).
+Assumes profile-mode data and a monotonic m/z axis.
+"""
+function detect_peaks_profile(mz::AbstractVector{<:Real},
+ y::AbstractVector{<:Real};
+ half_window::Int=10,
+ snr_threshold::Float64=2.0)
+ n = length(y)
+ n < 3 && return (Float64[], Float64[])
+
+ # Noise estimation using Median Absolute Deviation (robust to peaks)
+ noise_level = mad(y, normalize=true) + eps(Float64)
+
+ # Smooth the spectrum to make peak detection more robust
+ ys = smooth_spectrum(y; window=max(5, 2*half_window+1), order=2)
+
+ peak_idx = Int[]
+ @inbounds for i in 2:n-1
+ left = max(1, i - half_window)
+ right = min(n, i + half_window)
+
+ local_max = ys[i]
+ # A point is a peak if it's the maximum in its neighborhood and above the SNR threshold
+ if local_max >= maximum( @view ys[left:right]) && (local_max > snr_threshold * noise_level)
+ # Ensure we only record one point for flat-topped peaks
+ if isempty(peak_idx) || (i - last(peak_idx) > half_window)
+ push!(peak_idx, i)
+ end
+ end
+ end
+
+ # Return original intensities at peak locations
+ pk_mz = [float(mz[i]) for i in peak_idx]
+ pk_int = [float(y[i]) for i in peak_idx]
+
+ return (pk_mz, pk_int)
+end
+
+"""
+ detect_peaks_centroid(mz, y; intensity_threshold=0.0)
+
+Filters centroided data based on a minimum intensity threshold.
+"""
+function detect_peaks_centroid(mz::AbstractVector{<:Real},
+ y::AbstractVector{<:Real};
+ intensity_threshold::Float64=0.0)
+
+ keep_indices = findall(y .>= intensity_threshold)
+
+ return (mz[keep_indices], y[keep_indices])
+end
+
+# =============================================================================
+# 5) Peak Alignment
+# =============================================================================
+
+"""
+ align_peaks_lowess(ref_mz, tgt_mz; tolerance=0.002) -> warp::Function
+
+Generates a warping function `warp(x)` to map target m/z values to reference m/z values.
+Uses a lightweight LOWESS-like approach with linear interpolation.
+"""
+function align_peaks_lowess(ref_mz::Vector{<:Real},
+ tgt_mz::Vector{<:Real};
+ tolerance::Float64=0.002)
+ # Efficiently match peaks between sorted lists
+ pairs = Tuple{Float64,Float64}[] # (target_mz, reference_mz)
+ i = 1; j = 1
+ while i <= length(tgt_mz) && j <= length(ref_mz)
+ dt = tgt_mz[i] - ref_mz[j]
+ if abs(dt) <= tolerance
+ push!(pairs, (float(tgt_mz[i]), float(ref_mz[j])))
+ i += 1; j += 1
+ elseif dt < 0
+ i += 1
+ else
+ j += 1
+ end
+ end
+
+ if length(pairs) < 3
+ @warn "Too few matching peaks for alignment. Returning identity function."
+ return x -> float.(x)
+ end
+
+ t = [p[1] for p in pairs]
+ r = [p[2] for p in pairs]
+
+ # Lightly smooth the mapping to reduce noise
+ t_s = smooth_spectrum(t; window=5, order=2)
+ r_s = smooth_spectrum(r; window=5, order=2)
+
+ # Return a function that performs linear interpolation for warping
+ function warp(x::AbstractVector{<:Real})
+ out = similar(collect(float.(x)))
+ for (k, xv) in enumerate(x)
+ if xv <= t_s[1]
+ # Linear extrapolation at the start
+ m = (r_s[2]-r_s[1]) / (t_s[2]-t_s[1] + eps())
+ out[k] = r_s[1] + m*(xv - t_s[1])
+ elseif xv >= t_s[end]
+ # Linear extrapolation at the end
+ m = (r_s[end]-r_s[end-1]) / (t_s[end]-t_s[end-1] + eps())
+ out[k] = r_s[end-1] + m*(xv - t_s[end-1])
+ else
+ # Linear interpolation for points in the middle
+ lo = searchsortedlast(t_s, xv)
+ hi = lo + 1
+ α = (xv - t_s[lo]) / (t_s[hi] - t_s[lo] + eps())
+ out[k] = (1-α)*r_s[lo] + α*r_s[hi]
+ end
+ end
+ return out
+ end
+
+ return warp
+end
+
+# =============================================================================
+# 6) Peak Binning & Feature Matrix Generation
+# =============================================================================
+
+"""
+ _find_bin_index(x, bins) -> Int
+
+Efficiently finds the index of the bin `(low, high)` that contains `x` using binary search.
+Returns 0 if not found.
+"""
+function _find_bin_index(x::Float64, bins::Vector{Tuple{Float64,Float64}})
+ lo, hi = 1, length(bins)
+ while lo <= hi
+ mid = (lo + hi) >>> 1
+ b = bins[mid]
+ if x < b[1]
+ hi = mid - 1
+ elseif x > b[2]
+ lo = mid + 1
+ else
+ return mid
+ end
+ end
+ return 0
+end
+
+"""
+ bin_peaks(all_pk_mz, all_pk_int, tolerance; frequency_threshold=0.25)
+
+Groups peaks from all spectra into consensus m/z bins and creates a feature matrix.
+Filters out features that do not appear in a minimum fraction of spectra.
+"""
+function bin_peaks(all_pk_mz::Vector{<:AbstractVector{<:Real}},
+ all_pk_int::Vector{<:AbstractVector{<:Real}},
+ tolerance::Float64; frequency_threshold::Float64=0.25)
+ ns = length(all_pk_mz)
+ ns == 0 && return (zeros(0,0), Tuple{Float64,Float64}[])
+
+ # 1) Collect all unique peak m/z values and sort them
+ flat_mz = Float64[]
+ for v in all_pk_mz
+ append!(flat_mz, float.(v))
+ end
+ sort!(flat_mz)
+ isempty(flat_mz) && return (zeros(ns, 0), Tuple{Float64,Float64}[])
+
+ # 2) Create contiguous m/z bins based on tolerance
+ bins = Tuple{Float64,Float64}[]
+ cur_lo = flat_mz[1]
+ cur_hi = flat_mz[1]
+ for x in @view flat_mz[2:end]
+ if x - cur_hi <= tolerance
+ cur_hi = x # Extend the current bin
+ else
+ push!(bins, (cur_lo, cur_hi)) # Finalize old bin
+ cur_lo = x; cur_hi = x # Start a new one
+ end
+ end
+ push!(bins, (cur_lo, cur_hi))
+
+ # 3) Create the feature matrix (samples x features) using max intensity per bin
+ X = zeros(Float64, ns, length(bins))
+ for i in 1:ns
+ for (mzv, iv) in zip(all_pk_mz[i], all_pk_int[i])
+ bidx = _find_bin_index(float(mzv), bins)
+ if bidx > 0
+ X[i, bidx] = max(X[i, bidx], float(iv))
+ end
+ end
+ end
+
+ # 4) Filter features by minimum frequency
+ if frequency_threshold > 0
+ present_count = vec(sum(X .> 0, dims=1))
+ min_count = ceil(Int, frequency_threshold * ns)
+ keep_mask = findall(present_count .>= min_count)
+ X = X[:, keep_mask]
+ bins = bins[keep_mask]
+ end
+
+ return (X, bins)
+end
+
+# =============================================================================
+# 7) Plotting Helper
+# =============================================================================
+
+"""
+ plot_stage_spectrum(mz, intensity; title, ...)
+
+Returns a `CairoMakie.Figure` for a single spectrum trace. The caller is responsible for saving.
+"""
+function plot_stage_spectrum(mz::AbstractVector, intensity::AbstractVector;
+ title::AbstractString, xlabel::AbstractString="m/z",
+ ylabel::AbstractString="Intensity")
+ # This dynamic import is for script-like use; in a package, Makie would be a full dependency.
+ @eval begin
+ import CairoMakie
+ using CairoMakie
+ end
+ fig = CairoMakie.Figure(size = (1400, 500))
+ ax = CairoMakie.Axis(fig[1, 1], title=title, xlabel=xlabel, ylabel=ylabel)
+ CairoMakie.lines!(ax, mz, intensity)
+ return fig
+end
+
+# =============================================================================
+# 8) Pipeline Orchestrator
+# =============================================================================
+
+"""
+ run_preprocessing_pipeline(spectra; steps, params, on_stage)
+
+Executes a flexible preprocessing pipeline on a vector of spectra.
+
+# Arguments
+- `spectra`: A vector of `(mz, intensity)` tuples.
+- `steps`: A vector of symbols defining the pipeline order (e.g., `[:qc, :smooth, :baseline, :peaks, :bin]`).
+- `params`: A dictionary of parameters for each step.
+- `on_stage`: An optional callback function `on_stage(stage_symbol; idx, mz, intensity)` executed after each step for logging or visualization.
+
+# Returns
+- A `FeatureMatrix` if `:bin` is in the steps, otherwise the vector of processed spectra.
+"""
+function run_preprocessing_pipeline(spectra::Vector;
+ steps::Vector{Symbol},
+ params::Dict=Dict(),
+ on_stage::Function=(;kwargs...)->nothing)
+
+ processed = deepcopy(spectra) # Don't mutate the original input
+ reference_peaks = nothing # For alignment
+
+ _emit(stage::Symbol, idx::Int, mz, y) = on_stage(stage; idx=idx, mz=mz, intensity=y)
+
+ for step in steps
+ @info "Running step: $step"
+
+ if step === :qc
+ for (i, (mz, y)) in enumerate(processed)
+ (isempty(mz) || isempty(y)) && continue
+ _emit(:qc_raw, i, mz, y)
+ qc_is_empty(mz, y) && @warn "Spectrum at index $i is empty."
+ !qc_is_regular(mz) && @warn "m/z axis at index $i is not monotonic."
+ end
+
+ elseif step === :transform
+ meth = get(params, :transform_method, :sqrt)
+ for i in eachindex(processed)
+ mz, y = processed[i]
+ y_new = transform_intensity(y; method=meth)
+ processed[i] = (mz, y_new)
+ _emit(:transform, i, mz, y_new)
+ end
+
+ elseif step === :smooth
+ win = get(params, :sg_window, 21)
+ ord = get(params, :sg_order, 2)
+ for i in eachindex(processed)
+ mz, y = processed[i]
+ y_smooth = smooth_spectrum(y; window=win, order=ord)
+ processed[i] = (mz, y_smooth)
+ _emit(:smooth, i, mz, y_smooth)
+ end
+
+ elseif step === :baseline
+ iters = get(params, :snip_iterations, 100)
+ for i in eachindex(processed)
+ mz, y = processed[i]
+ baseline = snip_baseline(y, iters)
+ y_corrected = max.(0.0, y .- baseline)
+ processed[i] = (mz, y_corrected)
+ _emit(:baseline, i, mz, y_corrected)
+ end
+
+ elseif step === :normalize
+ mode = get(params, :normalize_method, :tic)
+ if mode === :tic
+ for i in eachindex(processed)
+ mz, y = processed[i]
+ y_norm = tic_normalize(y)
+ processed[i] = (mz, y_norm)
+ _emit(:normalize, i, mz, y_norm)
+ end
+ elseif mode === :pqn
+ # Note: PQN assumes spectra are on a common m/z grid.
+ matrix = hcat([float.(p[2]) for p in processed]...)
+ matrix_norm = pqn_normalize(matrix)
+ for i in eachindex(processed)
+ mz, _ = processed[i]
+ processed[i] = (mz, view(matrix_norm, :, i))
+ _emit(:normalize, i, mz, processed[i][2])
+ end
+ end
+
+ elseif step === :peaks
+ peak_results = Vector{Tuple{Vector{Float64},Vector{Float64}}}(undef, length(processed))
+ hw = get(params, :peak_half_window, 10)
+ snr = get(params, :peak_snr, 2.0)
+ for (i, (mz, y)) in enumerate(processed)
+ pk_mz, pk_int = detect_peaks_profile(mz, y; half_window=hw, snr_threshold=snr)
+ peak_results[i] = (pk_mz, pk_int)
+ # Emit with original mz axis but maybe stem plot of peaks?
+ _emit(:peaks, i, pk_mz, pk_int)
+ end
+ processed = peak_results
+ # Set reference for alignment
+ reference_peaks = isempty(processed) ? nothing : processed[1][1]
+
+ elseif step === :align
+ reference_peaks === nothing && (@error "Alignment requires a :peaks step first."; continue)
+ tol = get(params, :align_tolerance, 0.002)
+ for i in 2:length(processed)
+ tgt_peaks, intens = processed[i]
+ warp_func = align_peaks_lowess(reference_peaks, tgt_peaks; tolerance=tol)
+ processed[i] = (warp_func(tgt_peaks), intens)
+ _emit(:align, i, processed[i][1], processed[i][2])
+ end
+
+ elseif step === :bin
+ all_pks = [s[1] for s in processed]
+ all_ints = [s[2] for s in processed]
+ tol = get(params, :bin_tolerance, 0.002)
+ freq = get(params, :bin_min_frequency, 0.25)
+ mat, mz_bins = bin_peaks(all_pks, all_ints, tol; frequency_threshold=freq)
+ return FeatureMatrix(mat, mz_bins, collect(1:length(processed)))
+ end
+ end
+
+ @warn "Pipeline finished without a :bin step; returning processed spectra."
+ return processed
+end
+
+"""
+ run_preprocessing_pipeline(msi_data::MSIData, indices::Vector{Int}; steps, params, on_stage)
+
+Executes a flexible preprocessing pipeline on a subset of spectra from an MSIData object,
+with mode-aware logic for centroid and profile data.
+"""
+function run_preprocessing_pipeline(msi_data::MSIData, indices::Vector{Int};
+ steps::Vector{Symbol},
+ params::Dict=Dict(),
+ on_stage::Function=(;kwargs...)->nothing)
+
+ # This version of the pipeline is mode-aware.
+ # It processes spectra directly from the MSIData object.
+
+ processed_spectra = Vector{Tuple}(undef, length(indices))
+
+ # First, load all spectra and apply initial steps that run on individual spectra
+ for (i, spec_idx) in enumerate(indices)
+ mz, intensity = GetSpectrum(msi_data, spec_idx)
+ mode = msi_data.spectra_metadata[spec_idx].mode
+
+ on_stage(:qc_raw; idx=spec_idx, mz=mz, intensity=intensity)
+
+ for step in steps
+ if step === :transform
+ meth = get(params, :transform_method, :sqrt)
+ intensity = transform_intensity(intensity; method=meth)
+ on_stage(:transform; idx=spec_idx, mz=mz, intensity=intensity)
+ elseif step === :smooth && mode == PROFILE
+ win = get(params, :sg_window, 21)
+ ord = get(params, :sg_order, 2)
+ intensity = smooth_spectrum(intensity; window=win, order=ord)
+ on_stage(:smooth; idx=spec_idx, mz=mz, intensity=intensity)
+ elseif step === :baseline && mode == PROFILE
+ iters = get(params, :snip_iterations, 100)
+ baseline = snip_baseline(intensity, iters)
+ intensity = max.(0.0, intensity .- baseline)
+ on_stage(:baseline; idx=spec_idx, mz=mz, intensity=intensity)
+ elseif step === :normalize
+ norm_mode = get(params, :normalize_method, :tic)
+ if norm_mode === :tic
+ intensity = tic_normalize(intensity)
+ on_stage(:normalize; idx=spec_idx, mz=mz, intensity=intensity)
+ end
+ end
+ end
+ processed_spectra[i] = (mz, intensity, mode) # Store mode for peak detection
+ end
+
+ # Now, handle steps that require all spectra (like PQN) or are the final steps
+ final_result = nothing
+ for step in steps
+ if step === :normalize && get(params, :normalize_method, :tic) === :pqn
+ # Note: PQN assumes spectra are on a common m/z grid.
+ matrix = hcat([float.(p[2]) for p in processed_spectra]...)
+ matrix_norm = pqn_normalize(matrix)
+ for i in eachindex(processed_spectra)
+ mz, _, mode = processed_spectra[i]
+ processed_spectra[i] = (mz, view(matrix_norm, :, i), mode)
+ on_stage(:normalize; idx=indices[i], mz=mz, intensity=processed_spectra[i][2])
+ end
+ elseif step === :peaks
+ peak_results = Vector{Tuple{Vector{Float64},Vector{Float64}}}(undef, length(processed_spectra))
+ hw = get(params, :peak_half_window, 10)
+ snr = get(params, :peak_snr, 2.0)
+ intensity_thresh = get(params, :peak_intensity_threshold, 0.0)
+
+ for (i, (mz, y, mode)) in enumerate(processed_spectra)
+ spec_idx = indices[i]
+ if mode == PROFILE
+ pk_mz, pk_int = detect_peaks_profile(mz, y; half_window=hw, snr_threshold=snr)
+ else # CENTROID
+ pk_mz, pk_int = detect_peaks_centroid(mz, y; intensity_threshold=intensity_thresh)
+ end
+ peak_results[i] = (pk_mz, pk_int)
+ on_stage(:peaks; idx=spec_idx, mz=pk_mz, intensity=pk_int)
+ end
+ processed_spectra = peak_results # Now contains peak lists
+
+ elseif step === :align
+ # Alignment requires a reference peak list, typically from the first spectrum
+ reference_peaks = isempty(processed_spectra) ? nothing : processed_spectra[1][1]
+ if reference_peaks === nothing
+ @error "Alignment requires a :peaks step first."; continue
+ end
+ tol = get(params, :align_tolerance, 0.002)
+ for i in 2:length(processed_spectra)
+ tgt_peaks, intens = processed_spectra[i]
+ warp_func = align_peaks_lowess(reference_peaks, tgt_peaks; tolerance=tol)
+ processed_spectra[i] = (warp_func(tgt_peaks), intens)
+ on_stage(:align; idx=indices[i], mz=processed_spectra[i][1], intensity=processed_spectra[i][2])
+ end
+
+ elseif step === :bin
+ all_pks = [s[1] for s in processed_spectra]
+ all_ints = [s[2] for s in processed_spectra]
+ tol = get(params, :bin_tolerance, 0.002)
+ freq = get(params, :bin_min_frequency, 0.25)
+ mat, mz_bins = bin_peaks(all_pks, all_ints, tol; frequency_threshold=freq)
+ final_result = FeatureMatrix(mat, mz_bins, indices)
+ break # Binning is the last step
+ end
+ end
+
+ if final_result !== nothing
+ return final_result
+ else
+ @warn "Pipeline finished without a :bin step; returning processed spectra."
+ return processed_spectra
+ end
+end
+
+# =============================================================================
+# 9) Quality Control Metrics
+# =============================================================================
+
+"""
+ calculate_ppm_error(measured_mz::Float64, theoretical_mz::Float64) -> Float64
+
+Calculates mass accuracy in parts-per-million (PPM).
+
+# Formula
+PPM = 10⁶ × |measured_mz - theoretical_mz| / theoretical_mz
+"""
+function calculate_ppm_error(measured_mz::Real, theoretical_mz::Real)
+ if theoretical_mz == 0
+ return Inf
+ end
+ return 1e6 * abs(Float64(measured_mz) - Float64(theoretical_mz)) / Float64(theoretical_mz)
+end
+
+"""
+ calculate_ppm_error_bulk(measured_mz::Vector{Float64}, theoretical_mz::Vector{Float64}) -> Vector{Float64}
+
+Calculates PPM errors for multiple mass values.
+"""
+function calculate_ppm_error_bulk(measured_mz::Vector{Real}, theoretical_mz::Vector{Real})
+ return [calculate_ppm_error(m, t) for (m, t) in zip(measured_mz, theoretical_mz)]
+end
+
+"""
+ calculate_resolution_fwhm(mz::Float64, profile_mz::Vector{Float64},
+ profile_intensity::Vector{Float64}) -> Float64
+
+Calculates mass resolution using Full Width at Half Maximum (FWHM).
+
+# Formula
+Resolution = m / Δm, where Δm is FWHM
+
+# Arguments
+- `mz`: Peak centroid m/z
+- `profile_mz`: Full m/z array from profile data
+- `profile_intensity`: Full intensity array from profile data
+
+# Returns
+Resolution or NaN if cannot be calculated
+"""
+function calculate_resolution_fwhm(mz::Real, profile_mz::AbstractVector{<:Real},
+ profile_intensity::AbstractVector{<:Real})
+
+ # Find peak center index
+ peak_idx = argmin(abs.(profile_mz .- mz))
+ peak_height = Float64(profile_intensity[peak_idx])
+ half_max = peak_height / 2
+
+ # Find left half-maximum point (interpolate for accuracy)
+ left_idx = find_last_below(profile_intensity[1:peak_idx], half_max)
+ if left_idx == 0 || left_idx == length(profile_intensity[1:peak_idx])
+ return NaN
+ end
+
+ # Linear interpolation for left FWHM
+ x1, x2 = Float64(profile_mz[left_idx]), Float64(profile_mz[left_idx+1])
+ y1, y2 = Float64(profile_intensity[left_idx]), Float64(profile_intensity[left_idx+1])
+ left_fwhm = x1 + (x2 - x1) * (half_max - y1) / (y2 - y1)
+
+ # Find right half-maximum point
+ right_slice = profile_intensity[peak_idx:end]
+ right_offset = find_first_below(right_slice, half_max)
+ if right_offset == 0 || right_offset == length(right_slice)
+ return NaN
+ end
+
+ right_idx = peak_idx + right_offset - 1
+ x1, x2 = Float64(profile_mz[right_idx-1]), Float64(profile_mz[right_idx])
+ y1, y2 = Float64(profile_intensity[right_idx-1]), Float64(profile_intensity[right_idx])
+ right_fwhm = x1 + (x2 - x1) * (half_max - y1) / (y2 - y1)
+
+ fwhm = right_fwhm - left_fwhm
+ return fwhm > 0 ? Float64(mz) / fwhm : NaN
+end
+
+# Helper functions for FWHM calculation
+function find_last_below(v::AbstractVector{<:Real}, threshold::Real)
+ for i in length(v):-1:2
+ if v[i] >= threshold && v[i-1] < threshold
+ return i-1
+ end
+ end
+ return 0
+end
+
+function find_first_below(v::AbstractVector{<:Real}, threshold::Real)
+ for i in 1:(length(v)-1)
+ if v[i] >= threshold && v[i+1] < threshold
+ return i+1
+ end
+ end
+ return 0
+end
+
+"""
+ analyze_mass_accuracy(msi_data, reference_peaks; ppm_tolerance=20.0)
+
+Analyzes mass accuracy across the dataset using known reference peaks.
+
+# Arguments
+- `msi_data`: Your MSI dataset
+- `reference_peaks`: Dict of theoretical m/z values -> compound names
+- `ppm_tolerance`: Initial tolerance for peak matching
+
+# Returns
+Comprehensive mass accuracy report
+"""
+function analyze_mass_accuracy(msi_data, reference_peaks::Dict{Float64,String};
+ ppm_tolerance::Float64=5.0, sample_spectra=100)
+ println("\n[ MASS ACCURACY ANALYSIS ]")
+ println("PPM tolerance: ", ppm_tolerance, " ppm")
+ println("Spectra to sample: ", sample_spectra)
+
+ theoretical_mz = sort(collect(keys(reference_peaks)))
+ ppm_errors = Float64[]
+ matched_peaks = Tuple{Float64,Float64,String}[] # (theoretical, measured, compound)
+
+ # Sample spectra across the dataset
+ if sample_spectra >= length(msi_data.spectra_metadata)
+ spectrum_indices = 1:length(msi_data.spectra_metadata)
+ else
+ spectrum_indices = round.(Int, range(1, length(msi_data.spectra_metadata), length=sample_spectra))
+ end
+
+ for idx in spectrum_indices
+ mz, intensity = GetSpectrum(msi_data, idx)
+
+ # Detect peaks in this spectrum
+ detected_peaks, _ = detect_peaks_profile(mz, intensity, snr_threshold=3.0)
+
+ # Match detected peaks to reference peaks
+ for (i, theoretical) in enumerate(theoretical_mz)
+ # Find closest detected peak within tolerance
+ distances = abs.(detected_peaks .- theoretical)
+ if !isempty(distances)
+ min_idx = argmin(distances)
+ min_distance = distances[min_idx]
+
+ ppm_error = calculate_ppm_error(detected_peaks[min_idx], theoretical)
+
+ if ppm_error <= ppm_tolerance
+ push!(ppm_errors, ppm_error)
+ push!(matched_peaks, (theoretical, detected_peaks[min_idx], reference_peaks[theoretical]))
+ end
+ end
+ end
+ end
+
+ if isempty(ppm_errors)
+ @warn "No peaks matched within $ppm_tolerance ppm tolerance"
+ return (mean_ppm=NaN, std_ppm=NaN, min_ppm=NaN, max_ppm=NaN, optimal_ppm=NaN, n_matches=0, matched_peaks=[], all_ppm_errors=[])
+ end
+
+ # Calculate statistics
+ mean_ppm = mean(ppm_errors)
+ std_ppm = std(ppm_errors)
+ min_ppm = minimum(ppm_errors)
+ max_ppm = maximum(ppm_errors)
+
+ # Determine optimal ppm tolerance (mean + 3σ covers ~99.7% of peaks for normal distribution)
+ optimal_ppm = mean_ppm + 3 * std_ppm
+
+ return (
+ mean_ppm = mean_ppm,
+ std_ppm = std_ppm,
+ min_ppm = min_ppm,
+ max_ppm = max_ppm,
+ optimal_ppm = optimal_ppm,
+ n_matches = length(ppm_errors),
+ matched_peaks = matched_peaks,
+ all_ppm_errors = ppm_errors
+ )
+end
+
+"""
+ get_common_calibration_standards(standard_type::Symbol)
+
+Returns common calibration masses for different instrument types.
+
+# Supported standards
+- `:maldi_pos`: Common MALDI-TOF positive mode calibrants
+- `:maldi_neg`: Common MALDI-TOF negative mode calibrants
+- `:esi_pos`: ESI positive mode calibrants
+- `:lcms`: LC-MS commonly used standards
+"""
+function get_common_calibration_standards(standard_type::Symbol=:maldi_pos)
+ standards = Dict{Float64,String}()
+
+ if standard_type == :maldi_pos
+ standards = Dict(
+ 104.10754 => "C5H4N2 (Imidazole)",
+ 175.11995 => "C6H15O4P (Glycerophosphocholine fragment)",
+ 226.15687 => "C10H20NO4P (Phosphocholine)",
+ 322.04810 => "[Glu1]-Fibrinopeptide B fragment",
+ 379.09247 => "C12H22O11 (Sucrose)",
+ 515.32539 => "C26H52NO7P (PC(16:0/0:0))",
+ 622.02896 => "C20H12O5S2 (1-Hydroxypyrene-3,6,8-trisulfate)",
+ 757.39917 => "C37H74NO8P (PC(34:1))",
+ 1046.54198 => "Angiotensin I",
+ 1296.68477 => "ACTH clip 1-17",
+ 1570.67744 => "ACTH clip 18-39",
+ 2465.19829 => "ACTH clip 7-38"
+ )
+ elseif standard_type == :maldi_neg
+ standards = Dict(
+ 112.98563 => "C2F3O2 (Trifluoroacetate)",
+ 152.99568 => "C2F6S (Perfluoroethylsulfonate)",
+ 214.00166 => "C4F7O2 (Heptafluorobutyrate)",
+ 264.93278 => "C6F6 (Hexafluorobenzene)",
+ 362.96198 => "C8F15O2 (Perfluorooctanoate)",
+ 466.96714 => "C10F17O2S (Perfluorooctanesulfonate)"
+ )
+ elseif standard_type == :esi_pos
+ standards = Dict(
+ 118.08626 => "C5H12NO2 (Valine)",
+ 175.11900 => "C6H15O4P (Phosphocholine fragment)",
+ 524.26496 => "C23H48NO7P (LysoPC(16:0))",
+ 622.02896 => "C20H12O5S2 (Standard)",
+ 922.00980 => "C18H18O6N3S3 (Ultramark 1621)"
+ )
+ end
+
+ return standards
+end
+
+"""
+ generate_qc_report(msi_data; reference_peaks, output_dir)
+
+Generates a comprehensive QC report including mass accuracy and resolution.
+"""
+function generate_qc_report(msi_data, filename::String; reference_peaks=nothing, output_dir="qc_results", sample_spectra=100)
+ println("\n[ QC REPORT GENERATION ]")
+ println("Input file: ", filename)
+ println("Output directory: ", output_dir)
+ mkpath(output_dir)
+
+ # Use default calibrants if none provided
+ if reference_peaks === nothing
+ reference_peaks = get_common_calibration_standards(:maldi_pos)
+ end
+
+ println("Generating QC Report...")
+ println("Using $(length(reference_peaks)) reference masses")
+
+ # 1. Analyze mass accuracy
+ accuracy_report = analyze_mass_accuracy(msi_data, reference_peaks, sample_spectra=sample_spectra)
+
+ println("\n" * "="^50)
+ println("MASS ACCURACY REPORT")
+ println("="^50)
+ println("Mean PPM error: $(round(accuracy_report.mean_ppm, digits=2)) ppm")
+ println("Std PPM error: $(round(accuracy_report.std_ppm, digits=2)) ppm")
+ println("Min PPM error: $(round(accuracy_report.min_ppm, digits=2)) ppm")
+ println("Max PPM error: $(round(accuracy_report.max_ppm, digits=2)) ppm")
+ if haskey(accuracy_report, :optimal_ppm)
+ println("Optimal PPM tolerance: $(round(accuracy_report.optimal_ppm, digits=2)) ppm")
+ else
+ println("Optimal PPM tolerance: Not available")
+ end
+ println("Number of matches: $(accuracy_report.n_matches)")
+
+ # 2. Calculate resolution for a few representative peaks
+ println("\n" * "="^50)
+ println("RESOLUTION ANALYSIS")
+ println("="^50)
+
+ resolution_results = []
+ sample_spectra = min(10, length(msi_data.spectra_metadata))
+
+ for (i, idx) in enumerate(round.(Int, range(1, length(msi_data.spectra_metadata), length=sample_spectra)))
+ process_spectrum(msi_data, idx) do mz, intensity
+ if !qc_is_empty(mz, intensity)
+ # Test resolution on the most intense peak
+ max_intensity_idx = argmax(intensity)
+ test_mz = mz[max_intensity_idx]
+
+ resolution = calculate_resolution_fwhm(test_mz, mz, intensity)
+ if !isnan(resolution)
+ push!(resolution_results, resolution)
+ println("Spectrum $idx: Resolution = $(round(resolution))")
+ end
+ end
+ end
+ end
+
+ if !isempty(resolution_results)
+ avg_resolution = mean(resolution_results)
+ println("\nAverage resolution: $(round(avg_resolution))")
+ println("Resolution range: $(round(minimum(resolution_results))) - $(round(maximum(resolution_results)))")
+ end
+
+ # 3. Save detailed results
+
+ # Save PPM error distribution
+ ppm_df = DataFrame(
+ theoretical_mz = [p[1] for p in accuracy_report.matched_peaks],
+ measured_mz = [p[2] for p in accuracy_report.matched_peaks],
+ compound = [p[3] for p in accuracy_report.matched_peaks],
+ ppm_error = accuracy_report.all_ppm_errors
+ )
+
+ CSV.write(joinpath(output_dir, "mass_accuracy_results.csv"), ppm_df)
+
+ # Save resolution results
+ if !isempty(resolution_results)
+ res_df = DataFrame(resolution = resolution_results)
+ CSV.write(joinpath(output_dir, "resolution_results.csv"), res_df)
+ end
+
+ # 4. Create summary
+ summary = """
+ QC REPORT SUMMARY
+ =================
+ Date: $(now())
+ File: $(filename)
+ Spectra analyzed: $(length(msi_data.spectra_metadata))
+
+ MASS ACCURACY:
+ - Mean PPM: $(round(accuracy_report.mean_ppm, digits=2)) ppm
+ - Std PPM: $(round(accuracy_report.std_ppm, digits=2)) ppm
+ - Recommended tolerance: $(round(accuracy_report.optimal_ppm, digits=2)) ppm (mean + 3 * std)
+
+ RESOLUTION:
+ - Average: $(isempty(resolution_results) ? "N/A" : string(round(mean(resolution_results))))
+ - Range: $(isempty(resolution_results) ? "N/A" : "$(round(minimum(resolution_results))) - $(round(maximum(resolution_results)))")
+
+ RECOMMENDATIONS:
+ - Use $(round(accuracy_report.optimal_ppm, digits=2)) ppm for peak matching
+ - Instrument performance: $(accuracy_report.mean_ppm < 5 ? "Excellent" : accuracy_report.mean_ppm < 10 ? "Good" : "Needs calibration")
+ """
+
+ open(joinpath(output_dir, "qc_summary.txt"), "w") do f
+ write(f, summary)
+ end
+
+ println("\nQC report saved to: $output_dir")
+ return accuracy_report, resolution_results
+end
+
diff --git a/test/readme.md b/test/readme.md
index f80ee72..197af0d 100644
--- a/test/readme.md
+++ b/test/readme.md
@@ -37,21 +37,24 @@ Before running the tests, you must edit the `test/run_tests.jl` file to point to
cd /path/to/your/JuliaMSI
```
-2. **Execute the Test Script**:
+2. **Execute the Different Test Scripts**:
Run the following command from the project's root directory. This will install the necessary dependencies and run the tests.
```bash
julia --project=. test/run_tests.jl
```
+ ```bash
+ julia --project=. test/run_preprocessing.jl
+ ```
3. **Check the Results**:
The script will print its progress to the console. Any generated images (plots and image slices) will be saved in the `test/results/` directory.
## Test Case Configuration
-You can customize the test run by editing the variables in `test/run_tests.jl`.
+You can customize the test run by editing the variables in `test/run_tests.jl` or `test/run_preprocessing.jl` respectively.
### Enabling and Disabling Test Cases
-You can run or skip specific test cases by setting the corresponding boolean variables to `true` or `false`.
+You can run or skip specific test cases of run_tests by setting the corresponding boolean variables to `true` or `false`.
```julia
test1 = true # Runs Test Case 1
diff --git a/test/run_preprocessing.jl b/test/run_preprocessing.jl
new file mode 100644
index 0000000..4ea526a
--- /dev/null
+++ b/test/run_preprocessing.jl
@@ -0,0 +1,468 @@
+# test/run_preprocessing.jl
+
+# ===================================================================
+# Test Environment for the Preprocessing.jl Module
+# ===================================================================
+# This script tests the full preprocessing pipeline on single spectra
+# and total spectra from both .mzML and .imzML files.
+# It generates an overlay plot showing all preprocessing stages and
+# saves the resulting feature matrix to a CSV file.
+#
+# Instructions:
+# 1. Ensure the file paths in the "CONFIG" section are correct.
+# 2. Run the script from the project's root directory:
+# julia test/run_preprocessing.jl
+# 3. Check the `test/results/` folder for output plots and CSVs.
+# ===================================================================
+
+using Printf
+using CairoMakie
+import Pkg
+using DataFrames # For saving FeatureMatrix to CSV
+using CSV # For saving FeatureMatrix to CSV
+using Statistics # For mean()
+
+# --- Load Modules ---
+# Activate the project environment to access dependencies
+Pkg.activate(joinpath(@__DIR__, ".."))
+using MSI_src # This brings in Preprocessing.jl functions via export
+
+# ===================================================================
+# CONFIG: Test files and parameters
+# ===================================================================
+
+# --- Test Files ---
+# An mzML file for testing spectrum-based processing
+# const TEST_MZML_FILE = "/home/pixel/Documents/Cinvestav_2025/Analisis/CE4_BF_R1/CE4_BF_R1.mzML"
+# const TEST_MZML_FILE = "/home/pixel/Documents/Cinvestav_2025/Analisis/set de datos MS/Leaf_profile_LD_LTP_MS.mzML"
+const TEST_MZML_FILE = "/home/pixel/Documents/Cinvestav_2025/Analisis/set de datos MS/Escopolamina_tuneo_fraq_20ev.mzML"
+#const TEST_MZML_FILE = "/home/pixel/Documents/Cinvestav_2025/Analisis/set de datos MS/Atropina_tuneo_fraq_20ev.mzML"
+
+const MZML_SPECTRUM_ID = 1
+
+# An imzML file for testing
+# const TEST_IMZML_FILE = "/home/pixel/Documents/Cinvestav_2025/Analisis/CE4_BF_R1/CE4_BF_R1.imzML"
+# const IMZML_COORDS = (50, 50)
+const TEST_IMZML_FILE = "/home/pixel/Documents/Cinvestav_2025/Analisis/salida/Stomach_DHB_uncompressed.imzML"
+const IMZML_COORDS = (1997, 639)
+
+# --- Output Directory ---
+const RESULTS_DIR = "test/results"
+
+# ===================================================================
+# HELPER FUNCTIONS FOR PLOTTING
+# ===================================================================
+
+"""
+ plot_overlay_stages(collected_data, output_path, title)
+
+Creates a single plot overlaying spectra from different preprocessing stages.
+"""
+function plot_overlay_stages(collected_data, output_path, title)
+ fig = Figure(size = (1400, 800))
+ ax = Axis(fig[1, 1], title=title, xlabel="m/z", ylabel="Intensity")
+
+ colors = Makie.wong_colors() # A good set of distinct colors
+
+ for (i, (stage, mz, intensity)) in enumerate(collected_data)
+ color = colors[mod1(i, length(colors))] # Cycle through colors
+
+ # Plot the spectrum as a line
+ lines!(ax, mz, intensity, color=color, label=string(stage))
+
+ # If it's the peaks stage, also mark the peak tops
+ if stage == :peaks
+ scatter!(ax, mz, intensity, color=color, marker=:circle, markersize=8, label="$(string(stage)) (tops)")
+ end
+ end
+ axislegend(ax, position=:rt) # Right top position
+ save(output_path, fig)
+ println("SUCCESS: Overlay plot saved to $output_path")
+end
+
+
+# ===================================================================
+# TEST DEFINITIONS
+# ===================================================================
+
+"""
+ test_full_pipeline(msi_data, spectrum_id; output_dir, file_type_prefix)
+
+Tests the full preprocessing pipeline on a single spectrum and saves a plot
+for each intermediate step using the `on_stage` callback.
+`spectrum_id` can be an `Int` (for mzML) or a `Tuple{Int, Int}` (for imzML).
+"""
+function test_full_pipeline(msi_data, spectrum_id; output_dir, file_type_prefix, mz_tolerance=0.002)
+ println("\n--- Testing Full Preprocessing Pipeline on Spectrum: $spectrum_id (File Type: $file_type_prefix) ---")
+
+ # 1. Determine the spectrum index
+ local spec_idx
+ if spectrum_id isa Int
+ spec_idx = spectrum_id
+ else # Tuple for imzML
+ spec_idx = msi_data.coordinate_map[spectrum_id...]
+ end
+
+ if spec_idx == 0
+ println("SKIPPED: No spectrum found at coordinates $spectrum_id.")
+ return
+ end
+
+ # 2. Define the pipeline steps in the desired order
+ pipeline_steps = [
+ :qc,
+ :transform,
+ :smooth,
+ :baseline,
+ :normalize,
+ :peaks,
+ :align, # Align requires multiple spectra, but we'll run it on a single one for now (will warn)
+ :bin
+ ]
+
+ # Define parameters for each step
+ params = Dict(
+ :transform_method => :sqrt,
+ :sg_window => 15,
+ :sg_order => 2,
+ :snip_iterations => 100,
+ :normalize_method => :tic,
+ :peak_half_window => 10,
+ :peak_snr => 3.0,
+ :peak_intensity_threshold => 0.0, # For centroid peak detection
+ :align_tolerance => mz_tolerance,
+ :bin_tolerance => mz_tolerance,
+ :bin_min_frequency => 0.0 # Keep all bins for a single spectrum
+ )
+
+ # 3. Define the on_stage callback to collect data for overlay plot and save separate plots
+ collected_stage_data = []
+ stage_counter = Ref(0) # Initialize counter for sequential naming
+ normalized_spectrum = nothing # Variable to hold the normalized spectrum
+
+ function stage_callback(stage; idx, mz, intensity)
+ stage_counter[] += 1 # Increment counter
+ println(" -> Generating plot for stage: $stage")
+
+ local fig # Make fig available in the whole function scope
+
+ if stage == :normalize
+ normalized_spectrum = (mz, intensity)
+ fig = plot_stage_spectrum(mz, intensity, title="Stage: $stage (Spectrum $spectrum_id)")
+ elseif stage == :peaks && normalized_spectrum !== nothing
+ # For the peaks stage, plot the normalized spectrum as a base layer
+ fig = Figure(size = (1400, 500))
+ ax = Axis(fig[1, 1], title="Stage: Peaks (Spectrum $spectrum_id)", xlabel="m/z", ylabel="Intensity")
+ lines!(ax, normalized_spectrum[1], normalized_spectrum[2], color=:gray, label="Normalized Spectrum")
+ scatter!(ax, mz, intensity, color=:red, marker=:circle, markersize=8, label="Detected Peaks")
+ axislegend(ax)
+ else
+ # Default plotting for all other stages
+ fig = plot_stage_spectrum(mz, intensity, title="Stage: $stage (Spectrum $spectrum_id)")
+ end
+
+ # Save the figure
+ stage_output_path = joinpath(output_dir, "$(file_type_prefix)_$(spectrum_id)_$(stage_counter[])_$(stage).png")
+ save(stage_output_path, fig)
+
+ # Collect data for overlay plot
+ push!(collected_stage_data, (stage, mz, intensity))
+ end
+
+ # 4. Run the pipeline on the single spectrum
+ println("Running pipeline with steps: $pipeline_steps")
+ processed_result = run_preprocessing_pipeline(
+ msi_data,
+ [spec_idx], # The pipeline expects a vector of indices
+ steps=pipeline_steps,
+ params=params,
+ on_stage=stage_callback
+ )
+
+ # 5. Generate and save the overlay plot
+ overlay_output_path = joinpath(output_dir, "$(file_type_prefix)_$(spectrum_id)_all_stages_overlay.png")
+ plot_overlay_stages(collected_stage_data, overlay_output_path, "Preprocessing Stages Overlay (Spectrum $spectrum_id)")
+
+ # 6. Save feature matrix if generated
+ if processed_result isa FeatureMatrix
+ feature_matrix_output_path = joinpath(output_dir, "$(file_type_prefix)_$(spectrum_id)_feature_matrix.csv")
+ # Convert mz_bins to a more readable format for CSV
+ mz_labels = ["$(round(b[1], digits=4))_$(round(b[2], digits=4))" for b in processed_result.mz_bins]
+ df = DataFrame(processed_result.matrix, Symbol.(mz_labels))
+ CSV.write(feature_matrix_output_path, df)
+ println("SUCCESS: Feature matrix saved to $feature_matrix_output_path")
+ else
+ @warn "Pipeline did not return a FeatureMatrix for Spectrum $spectrum_id."
+ processed_result
+ end
+
+ println("--- Pipeline test finished for Spectrum: $spectrum_id (File Type: $file_type_prefix) ---")
+ println("Check the '$(output_dir)' directory for output plots and CSVs.")
+end
+
+"""
+ test_full_pipeline_on_total_spectrum(msi_data; output_dir, file_type_prefix)
+
+Tests the full preprocessing pipeline on the *total spectrum* (sum of all spectra)
+and saves a plot for each intermediate step.
+"""
+function test_full_pipeline_on_total_spectrum(msi_data; output_dir, file_type_prefix, mz_tolerance=0.002)
+ println("\n--- Testing Full Preprocessing Pipeline on TOTAL Spectrum (File Type: $file_type_prefix) ---")
+
+ # 1. Get the total spectrum
+ total_mz, total_intensity = get_total_spectrum(msi_data)
+ total_spectrum = (total_mz, total_intensity)
+
+ if qc_is_empty(total_mz, total_intensity)
+ println("SKIPPED: Total spectrum is empty.")
+ return
+ end
+
+ # 2. Define the pipeline steps and parameters (same as for single spectrum)
+ pipeline_steps = [
+ :qc,
+ :transform,
+ :smooth,
+ :baseline,
+ :normalize,
+ :peaks,
+ :align, # Align requires multiple spectra, but we'll run it on a single one for now (will warn)
+ :bin
+ ]
+
+ params = Dict(
+ :transform_method => :sqrt,
+ :sg_window => 15,
+ :sg_order => 2,
+ :snip_iterations => 100,
+ :normalize_method => :tic,
+ :peak_half_window => 10,
+ :peak_snr => 3.0,
+ :align_tolerance => mz_tolerance,
+ :bin_tolerance => mz_tolerance,
+ :bin_min_frequency => 0.0 # Keep all bins for a single spectrum
+ )
+
+ # 3. Define the on_stage callback
+ collected_stage_data = []
+ stage_counter = Ref(0) # Initialize counter for sequential naming
+ normalized_spectrum_total = nothing # Variable to hold the normalized spectrum
+
+ function stage_callback_total(stage; idx, mz, intensity)
+ stage_counter[] += 1 # Increment counter
+ println(" -> Generating plot for stage: $stage (Total Spectrum)")
+
+ local fig
+
+ if stage == :normalize
+ normalized_spectrum_total = (mz, intensity)
+ fig = plot_stage_spectrum(mz, intensity, title="Stage: $stage (Total Spectrum)")
+ elseif stage == :peaks && normalized_spectrum_total !== nothing
+ fig = Figure(size = (1400, 500))
+ ax = Axis(fig[1, 1], title="Stage: Peaks (Total Spectrum)", xlabel="m/z", ylabel="Intensity")
+ lines!(ax, normalized_spectrum_total[1], normalized_spectrum_total[2], color=:gray, label="Normalized Spectrum")
+ scatter!(ax, mz, intensity, color=:red, marker=:circle, markersize=8, label="Detected Peaks")
+ axislegend(ax)
+ else
+ fig = plot_stage_spectrum(mz, intensity, title="Stage: $stage (Total Spectrum)")
+ end
+
+ # Save separate plot with sequential name
+ stage_output_path = joinpath(output_dir, "$(file_type_prefix)_total_$(stage_counter[])_$(stage).png")
+ save(stage_output_path, fig)
+
+ # Collect data for overlay plot
+ push!(collected_stage_data, (stage, mz, intensity))
+ end
+
+ # 4. Run the pipeline on the single total spectrum
+ println("Running pipeline with steps: $pipeline_steps")
+ processed_result = run_preprocessing_pipeline(
+ [total_spectrum], # Pass the total spectrum as a vector of one spectrum
+ steps=pipeline_steps,
+ params=params,
+ on_stage=stage_callback_total
+ )
+
+ # 5. Generate and save the overlay plot
+ overlay_output_path = joinpath(output_dir, "$(file_type_prefix)_total_all_stages_overlay.png")
+ plot_overlay_stages(collected_stage_data, overlay_output_path, "Preprocessing Stages Overlay (Total Spectrum)")
+
+ # 6. Save feature matrix if generated
+ if processed_result isa FeatureMatrix
+ feature_matrix_output_path = joinpath(output_dir, "$(file_type_prefix)_total_feature_matrix.csv")
+ # Convert mz_bins to a more readable format for CSV
+ mz_labels = ["$(round(b[1], digits=4))_$(round(b[2], digits=4))" for b in processed_result.mz_bins]
+ df = DataFrame(processed_result.matrix, Symbol.(mz_labels))
+ CSV.write(feature_matrix_output_path, df)
+ println("SUCCESS: Feature matrix saved to $feature_matrix_output_path")
+ else
+ @warn "Pipeline did not return a FeatureMatrix for Total Spectrum."
+ processed_result
+ end
+
+ println("--- Pipeline test finished for TOTAL Spectrum (File Type: $file_type_prefix) ---")
+ println("Check the '$(output_dir)' directory for output plots and CSVs.")
+end
+
+
+# ===================================================================
+# TEST RUNNER
+# ===================================================================
+
+function run_preprocessing_tests()
+ println("="^80)
+ println("STARTING PREPROCESSING TEST SUITE")
+ println("="^80)
+
+ # --- Test Case 1: Run full pipeline on a single mzML spectrum ---
+ println("\n" * "="^20 * " Test Case 1: Full Pipeline on .mzML Spectrum " * "="^20)
+ println("FILE: ", TEST_MZML_FILE)
+ if isfile(TEST_MZML_FILE)
+ try
+ msi_data_mzml = OpenMSIData(TEST_MZML_FILE)
+
+ # Dynamically determine tolerance
+ println("\n--- Calculating optimal tolerance for .mzML data ---")
+ report_mzml = analyze_mass_accuracy(msi_data_mzml, get_common_calibration_standards(:maldi_pos))
+ mz_tolerance_mzml = 0.002 # Default
+ if haskey(report_mzml, :optimal_ppm) && !isnan(report_mzml.optimal_ppm) && !isempty(report_mzml.matched_peaks)
+ avg_mz = mean([p[1] for p in report_mzml.matched_peaks])
+ mz_tolerance_mzml = avg_mz * report_mzml.optimal_ppm / 1e6
+ println("Optimal PPM: $(round(report_mzml.optimal_ppm, digits=2)), Average m/z: $(round(avg_mz, digits=2))")
+ println("Calculated m/z tolerance: $(round(mz_tolerance_mzml, digits=5))")
+ else
+ println("Could not determine optimal tolerance, using default: $mz_tolerance_mzml")
+ end
+
+ # Create a dedicated subdirectory for the output plots
+ mzml_output_dir = joinpath(RESULTS_DIR, "mzml_pipeline_stages")
+ mkpath(mzml_output_dir)
+
+ test_full_pipeline(msi_data_mzml, MZML_SPECTRUM_ID, output_dir=mzml_output_dir, file_type_prefix="mzml", mz_tolerance=mz_tolerance_mzml)
+ test_full_pipeline_on_total_spectrum(msi_data_mzml, output_dir=mzml_output_dir, file_type_prefix="mzml", mz_tolerance=mz_tolerance_mzml)
+ catch e
+ println("ERROR in .mzML pipeline test: $e")
+ showerror(stdout, e, catch_backtrace())
+ end
+ else
+ println("SKIPPED: File not found: $TEST_MZML_FILE")
+ end
+
+ # --- Test Case 2: Run full pipeline on a single imzML spectrum ---
+ println("\n" * "="^20 * " Test Case 2: Full Pipeline on .imzML Spectrum " * "="^20)
+ println("FILE: ", TEST_IMZML_FILE)
+ if isfile(TEST_IMZML_FILE)
+ try
+ msi_data_imzml = OpenMSIData(TEST_IMZML_FILE)
+
+ # Dynamically determine tolerance
+ println("\n--- Calculating optimal tolerance for .imzML data ---")
+ report_imzml = analyze_mass_accuracy(msi_data_imzml, get_common_calibration_standards(:maldi_pos))
+ mz_tolerance_imzml = 0.002 # Default
+ if haskey(report_imzml, :optimal_ppm) && !isnan(report_imzml.optimal_ppm) && !isempty(report_imzml.matched_peaks)
+ avg_mz = mean([p[1] for p in report_imzml.matched_peaks])
+ mz_tolerance_imzml = avg_mz * report_imzml.optimal_ppm / 1e6
+ println("Optimal PPM: $(round(report_imzml.optimal_ppm, digits=2)), Average m/z: $(round(avg_mz, digits=2))")
+ println("Calculated m/z tolerance: $(round(mz_tolerance_imzml, digits=5))")
+ else
+ println("Could not determine optimal tolerance, using default: $mz_tolerance_imzml")
+ end
+
+ # Create a dedicated subdirectory for the output plots
+ imzml_output_dir = joinpath(RESULTS_DIR, "imzml_pipeline_stages")
+ mkpath(imzml_output_dir)
+
+ test_full_pipeline(msi_data_imzml, IMZML_COORDS, output_dir=imzml_output_dir, file_type_prefix="imzml", mz_tolerance=mz_tolerance_imzml)
+ test_full_pipeline_on_total_spectrum(msi_data_imzml, output_dir=imzml_output_dir, file_type_prefix="imzml", mz_tolerance=mz_tolerance_imzml)
+ # generate_qc_report(msi_data_imzml, TEST_IMZML_FILE, output_dir=imzml_output_dir)
+ custom_reference_peaks = Dict(
+ 31.974 => "Red Phosphorus",
+ 432.6584 => "P13",
+ 464.6059 => "P15",
+ 526.5534 => "P17",
+ 650.4485 => "P21",
+ 774.3435 => "P25",
+ 898.2385 => "P29",
+ 950.1861 => "P31",
+ 1022.1336 => "P33",
+ 1146.0286 => "P37",
+ 1593.8187 => "P45",
+ 772.433 => "Unknown 1",
+ 772.5253 => "Unknown 2"
+ )
+ n_samples = length(msi_data_imzml.spectra_metadata)
+ generate_qc_report(msi_data_imzml, TEST_IMZML_FILE, reference_peaks=custom_reference_peaks, output_dir=imzml_output_dir, sample_spectra=n_samples)
+ catch e
+ println("ERROR in .imzML pipeline test: $e")
+ showerror(stdout, e, catch_backtrace())
+ end
+ else
+ println("SKIPPED: File not found: $TEST_IMZML_FILE")
+ end
+
+ # --- Test Case 3: Generate QC Report for .imzML data ---
+ println("\n" * "="^20 * " Test Case 3: QC Report Generation for .imzML " * "="^20)
+ println("FILE: ", TEST_IMZML_FILE)
+ if isfile(TEST_IMZML_FILE)
+ try
+ msi_data_imzml = OpenMSIData(TEST_IMZML_FILE)
+
+ # Create a dedicated subdirectory for the QC report
+ qc_output_dir = joinpath(RESULTS_DIR, "qc_report")
+ mkpath(qc_output_dir)
+
+ println("\n--- Generating comprehensive QC report ---")
+ custom_reference_peaks = Dict(
+ 31.974 => "Red Phosphorus",
+ 432.6584 => "P13",
+ 464.6059 => "P15",
+ 526.5534 => "P17",
+ 650.4485 => "P21",
+ 774.3435 => "P25",
+ 898.2385 => "P29",
+ 950.1861 => "P31",
+ 1022.1336 => "P33",
+ 1146.0286 => "P37",
+ 1593.8187 => "P45",
+ 772.433 => "Unknown 1",
+ 772.5253 => "Unknown 2"
+ )
+ # You can control the number of spectra sampled for the QC report.
+ # For the most accurate results, you can sample all spectra, but it will take longer.
+ # To sample all, use: n_samples = length(msi_data_imzml.spectra_metadata)
+ n_samples = length(msi_data_imzml.spectra_metadata)
+
+ #generate_qc_report(msi_data_imzml, TEST_IMZML_FILE, output_dir=qc_output_dir)
+ generate_qc_report(msi_data_imzml, TEST_IMZML_FILE, reference_peaks=custom_reference_peaks, output_dir=qc_output_dir, sample_spectra=n_samples)
+
+ println("\n--- Analyzing specific reference peaks ---")
+ #=
+ reference_peaks = Dict(
+ 104.10754 => "Imidazole",
+ 175.11995 => "GPC fragment",
+ 226.15687 => "Phosphocholine"
+ )
+ =#
+ report = analyze_mass_accuracy(msi_data_imzml, custom_reference_peaks)
+ if haskey(report, :optimal_ppm)
+ println("Optimal PPM tolerance with specific peaks: $(round(report.optimal_ppm, digits=2)) ppm")
+ else
+ println("Could not determine optimal PPM with specific peaks.")
+ end
+
+ catch e
+ println("ERROR in QC report generation test: $e")
+ showerror(stdout, e, catch_backtrace())
+ end
+ else
+ println("SKIPPED: File not found: $TEST_IMZML_FILE")
+ end
+
+ println("\nPreprocessing tests finished.")
+end
+
+# --- Execute ---
+# Ensure the results directory exists
+mkpath(RESULTS_DIR)
+@time run_preprocessing_tests()
\ No newline at end of file
diff --git a/test/run_tests.jl b/test/run_tests.jl
index c079ed9..212dc24 100644
--- a/test/run_tests.jl
+++ b/test/run_tests.jl
@@ -10,7 +10,7 @@
# Instructions:
# 1. Fill in the placeholder paths in the "CONFIG" section below.
# 2. Run the script from the project's root directory:
-# julia test/run_tests.jl
+# julia --project=. test/run_tests.jl
# 3. Check the `test/results/` folder for the output images.
# ===================================================================
+