diff --git a/Manifest.toml b/Manifest.toml index ff5e2fb..e1ae255 100644 --- a/Manifest.toml +++ b/Manifest.toml @@ -2,7 +2,7 @@ julia_version = "1.11.7" manifest_format = "2.0" -project_hash = "896c16630591d2732da8e3315c4b64c618fd5ebb" +project_hash = "0f3959248d174e05d23b9e9267b00c6d2ead256d" [[deps.ATK_jll]] deps = ["Artifacts", "Glib_jll", "JLLWrappers", "Libdl"] diff --git a/Project.toml b/Project.toml index 653695f..0db27f1 100644 --- a/Project.toml +++ b/Project.toml @@ -19,6 +19,8 @@ NativeFileDialog = "e1fe445b-aa65-4df4-81c1-2041507f0fd4" NaturalSort = "c020b1a1-e9b0-503a-9c33-f039bfc54a85" PlotlyBase = "a03496cd-edff-5a9b-9e67-9cda94a718b5" Printf = "de0858da-6303-5e67-8744-51eddeeeb8d7" +Serialization = "9e88b42a-f829-5b0c-bbe9-9e923198166b" Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2" +StatsBase = "2913bbd2-ae8a-5f71-8c99-4fb6c76f3a91" StipplePlotly = "ec984513-233d-481d-95b0-a3b58b97af2b" julia_mzML_imzML = "38eb50d3-2fb6-4afa-992a-964ed8562ed9" diff --git a/app.jl b/app.jl index b0fee09..b836bce 100644 --- a/app.jl +++ b/app.jl @@ -6,7 +6,8 @@ using Libz using PlotlyBase using CairoMakie using Colors -using julia_mzML_imzML +# using julia_mzML_imzML +using MSI_src # Import the new MSIData library using Statistics using NaturalSort using Images @@ -14,6 +15,10 @@ using LinearAlgebra using NativeFileDialog # Opens the file explorer depending on the OS using StipplePlotly using Base.Filesystem: mv # To rename files in the system + +# Bring MSIData into App module's scope +using .MSI_src: MSIData, OpenMSIData, GetSpectrum, IterateSpectra, ImzMLSource, _iterate_spectra_fast, MzMLSource, find_mass, ViridisPalette, get_mz_slice, quantize_intensity, save_bitmap, median_filter, save_bitmap, downsample_spectrum, TrIQ + include("./julia_imzML_visual.jl") @genietools @@ -57,6 +62,7 @@ include("./julia_imzML_visual.jl") @in compareBtn=false # To open dialog @in createMeanPlot=false # To generate mean spectrum plot @in createXYPlot=false # To generate an spectrum plot according to the xy values inputed + @in createSumPlot=false # To generate a sum of all the spectrum plots @in image3dPlot=false # To generate 3d plot based on current image @in triq3dPlot=false # To generate 3d plot based on current triq image @in imageCPlot=false # To generate contour plots of current image @@ -110,10 +116,12 @@ include("./julia_imzML_visual.jl") @out msgimgComp="" @out msgtriqComp="" + # Centralized MSIData object + @out msi_data::Union{MSIData, Nothing} = nothing + # Saves the route where imzML and mzML files are located @out full_route="" - @out full_routeMz="" - @out full_routeMz2="" + # For the creation of images with a more specific mass charge @out text_nmass="" @@ -145,10 +153,12 @@ include("./julia_imzML_visual.jl") layoutImg=PlotlyBase.Layout( xaxis=PlotlyBase.attr( visible=false, - scaleanchor="y" + scaleanchor="y", + range=[0, 0] ), yaxis=PlotlyBase.attr( - visible=false + visible=false, + range=[0, 0] ), margin=attr(l=0,r=0,t=0,b=0,pad=0) ) @@ -168,10 +178,10 @@ include("./julia_imzML_visual.jl") # Interface Plot Spectrum layoutSpectra=PlotlyBase.Layout( title="Spectrum plot", + hovermode="closest", xaxis=PlotlyBase.attr( title="m/z", - showgrid=true, - tickformat = ".3g" + showgrid=true ), yaxis=PlotlyBase.attr( title="Intensity", @@ -249,319 +259,305 @@ include("./julia_imzML_visual.jl") # Reactive handlers watch a variable and execute a block of code when its value changes # The onbutton handler will set the variable to false after the block is executed - @onbutton btnSearch begin - full_route=pick_file(; filterlist="imzML,imzml,mzML,mzml") - msg="" - if full_route != "" - if endswith(full_route, "imzml") - alt_route=replace(full_route, r"\.[^.]*$" => ".imzML") - mv(full_route, alt_route) - full_route=alt_route - # to detect if there's an mzml named wrong - alt_routeMz=replace(full_route, r"\.[^.]*$" => ".mzml") - if isfile(alt_routeMz) - alt_routeMz2=replace(alt_routeMz, r"\.[^.]*$" => ".mzML") - mv(alt_routeMz, alt_routeMz2) - alt_routeMz=alt_routeMz2 - end - end - if endswith(full_route, "mzml") - alt_route=replace(full_route, r"\.[^.]*$" => ".mzML") - mv(full_route, alt_route) - full_route = alt_route - # to detect if there's an imzml named wrong - alt_routeMz=replace(full_route, r"\.[^.]*$" => ".imzml") - if isfile(alt_routeMz) - alt_routeMz2=replace(alt_routeMz, r"\.[^.]*$" => ".imzML") - mv(alt_routeMz, alt_routeMz2) - alt_routeMz=alt_routeMz2 - end - end + @onbutton btnSearch @time begin + # This part is synchronous and blocking, which is unavoidable + picked_route = pick_file(; filterlist="imzML,imzml,mzML,mzml") + + if isempty(picked_route) + msg = "No file selected." + warning_msg = true + return end - if full_route=="" - msg="No file selected" - warning_msg=true - btnStartDisable=true - btnSpectraDisable=true - SpectraEnabled=false - else - if endswith(full_route, "imzML") # Case if the file loaded is imzML - btnStartDisable=false - btnPlotDisable=false - # Splitting the route with regex from imzml to mzml so the plotting can work - full_routeMz=replace(full_route, r"\.[^.]*$" => ".mzML") - if isfile(full_routeMz) - # Start the mean spectrum creation on loading - progressSpectraPlot=true - btnPlotDisable=true - btnStartDisable=true - msg="Loading mean spectrum plot..." - sTime=time() - plotdata, plotlayout, xSpectraMz, ySpectraMz=meanSpectrumPlot(full_routeMz) - selectedTab="tab2" - progressSpectraPlot=false - btnPlotDisable=false - # We enable coord search and spectra plot creation, also re-enable main process - btnSpectraDisable=false - SpectraEnabled=true - btnStartDisable=false - - fTime=time() - eTime=round(fTime-sTime,digits=3) - msg="Plot loaded in $(eTime) seconds" + # UI updates immediately + progress = true + msg = "Opening file: $(basename(picked_route))..." + + @async begin + try + # --- Normalize file extension and path --- + if endswith(picked_route, "imzml") + full_route = replace(picked_route, r"\.imzml$"i => ".imzML") + mv(picked_route, full_route, force=true) + elseif endswith(picked_route, "mzml") + full_route = replace(picked_route, r"\.mzml$"i => ".mzML") + mv(picked_route, full_route, force=true) else - # If there's no mzML file, we deny access again - btnSpectraDisable=true - SpectraEnabled=false + full_route = picked_route end - else # Case if the file loaded is mzML - full_routeMz=full_route - btnSpectraDisable=false - SpectraEnabled=true - if isfile(full_routeMz) - # Start the mean spectrum creation on loading - progressSpectraPlot=true - btnPlotDisable=true - btnStartDisable=true - msg="Loading mean spectrum plot..." - sTime=time() - plotdata, plotlayout, xSpectraMz, ySpectraMz=meanSpectrumPlot(full_routeMz) - selectedTab="tab2" - progressSpectraPlot=false - btnPlotDisable=false - # We enable coord search and spectra plot creation - btnSpectraDisable=false - SpectraEnabled=true + # --- Load data using the new MSIData library --- + sTime = time() + msi_data = OpenMSIData(full_route) + + w, h = msi_data.image_dims + imgWidth, imgHeight = w > 0 ? (w, h) : (500, 500) - fTime=time() - eTime=round(fTime-sTime,digits=3) - msg="Plot loaded in $(eTime) seconds" - end - # Splitting the route the same way - full_route=replace(full_route, r"\.[^.]*$" => ".imzML") - if isfile(full_route) - btnStartDisable=false - else - btnStartDisable=true - full_route=full_routeMz - end + fTime = time() + eTime = round(fTime - sTime, digits=3) + msg = "File loaded in $(eTime) seconds. Calculating total spectrum..." + + # Enable UI controls + btnStartDisable = !(msi_data.source isa ImzMLSource) + btnPlotDisable = false + btnSpectraDisable = false + SpectraEnabled = true + + # --- Automatically generate and display the sum spectrum plot --- + # This is still part of the same async task + progressSpectraPlot = true + sTime = time() + + plotdata, plotlayout, xSpectraMz, ySpectraMz = sumSpectrumPlot(msi_data) + + selectedTab = "tab2" + fTime = time() + eTime = round(fTime - sTime, digits=3) + msg = "Total spectrum plot loaded in $(eTime) seconds." + + catch e + msi_data = nothing + msg = "Error loading file: $e" + warning_msg = true + btnStartDisable = true + btnSpectraDisable = true + SpectraEnabled = false + @error "File loading failed" exception=(e, catch_backtrace()) + finally + # This now correctly runs after everything is finished + progress = false + progressSpectraPlot = false end - xCoord=0 - yCoord=0 end end - @onbutton mainProcess begin - progress=true # Start progress button animation - btnStartDisable=true # We disable the button to avoid multiple requests - btnPlotDisable=true - btnSpectraDisable=true - text_nmass=replace(string(Nmass), "." => "_") - sTime=time() - if isfile(full_route) && Nmass > 0 && Tol > 0 && Tol <=1 && colorLevel > 1 && colorLevel < 257 - msg="File exists, Nmass=$(Nmass) Tol=$(Tol). Loading file will begin, please be patient." + @onbutton mainProcess @time begin + # UI updates immediately + progress = true + btnStartDisable = true + btnPlotDisable = true + btnSpectraDisable = true + + @async begin try - spectra=LoadImzml(full_route) - msg="File loaded. Creating spectra with the specific mass and tolerance, please be patient." - slice=GetMzSliceJl(spectra,Nmass,Tol) - fig=CairoMakie.Figure(size=(150, 250)) # Container - # Append a query string to force the image to refresh - timestamp=string(time_ns()) - if triqEnabled # If we have TrIQ - if triqProb < 0.8 || triqProb > 1 - msg="Incorrect TrIQ values, please adjust accordingly and try again." - warning_msg=true - else - image_path=joinpath("./public", "TrIQ_$(text_nmass).bmp") - valid_slice=false - while Tol <= 1.0 && !valid_slice - try - slice=GetMzSliceJl(spectra, Nmass, Tol) - sliceTriq=TrIQ(slice, colorLevel, triqProb) - if MFilterEnabled # If the Median filter is ON - sliceTriq=medianFilterjl(sliceTriq) - end - valid_slice=true - catch e - msg="Warning: insufficient tolerance, inputs modified to allow the creation of an image regardless=$Tol: $e" - Tol += 0.1 - end - end - sliceTriq=reverse(sliceTriq, dims=2) - SaveBitmapCl(joinpath("public", "TrIQ_$(text_nmass).bmp"),sliceTriq,ViridisPalette) - # Use timestamp to refresh image interface container - imgIntT="/TrIQ_$(text_nmass).bmp?t=$(timestamp)" - plotdataImgT, plotlayoutImgT, imgWidth, imgHeight=loadImgPlot(imgIntT) - # Get current image - current_triq="TrIQ_$(text_nmass).bmp" - msgtriq="TrIQ image with the Nmass of $(replace(text_nmass, "_" => "."))" - # Create colorbar - bound=julia_mzML_imzML.GetOutlierThres(slice, triqProb) - levels=range(bound[1],stop=bound[2], length=8) - levels=vcat(levels, 2*levels[end]-levels[end-1]) - Colorbar(fig[1, 1], colormap=cgrad(:viridis, colorLevel, categorical=true), limits=(0, bound[2]),ticks=levels,tickformat=log_tick_formatter, label="Intensity", size=25) - save("public/colorbar_TrIQ_$(text_nmass).png", fig) - colorbarT="/colorbar_TrIQ_$(text_nmass).png?t=$(timestamp)" - # Get current colorbar - current_col_triq="colorbar_TrIQ_$(text_nmass).png" - # We update the directory to include the new placed images. - triq_bmp=sort(filter(filename -> startswith(filename, "TrIQ_") && endswith(filename, ".bmp"), readdir("public")),lt=natural) - col_triq_png=sort(filter(filename -> startswith(filename, "colorbar_TrIQ_") && endswith(filename, ".png"), readdir("public")),lt=natural) - fTime=time() - eTime=round(fTime-sTime,digits=3) - msg="The file has been created in $(eTime) seconds successfully inside the 'public' folder of the app" - selectedTab="tab1" - end - else # If we don't use TrIQ - image_path=joinpath("./public", "MSI_$(text_nmass).bmp") + text_nmass = replace(string(Nmass), "." => "_") + sTime = time() + + if msi_data === nothing || !(msi_data.source isa ImzMLSource) + msg = "No .imzML file loaded or selected file is not an .imzML. Please select a valid file." + warning_msg = true + elseif Nmass > 0 && Tol > 0 && Tol <= 1 && colorLevel > 1 && colorLevel < 257 + msg = "Creating image for Nmass=$(Nmass) Tol=$(Tol). Please be patient." try - sliceQuant=IntQuantCl(slice,Int(colorLevel-1)) - if MFilterEnabled # If the Median filter is ON - sliceQuant=medianFilterjl(sliceQuant) + # Use the new get_mz_slice with the centralized MSIData object + slice = get_mz_slice(msi_data, Nmass, Tol) + fig = CairoMakie.Figure(size=(150, 250)) # Container + timestamp = string(time_ns()) + + if triqEnabled # If we have TrIQ + if triqProb < 0.8 || triqProb > 1 + msg = "Incorrect TrIQ values, please adjust accordingly and try again." + warning_msg = true + else + sliceTriq = TrIQ(slice, colorLevel, triqProb) + if MFilterEnabled + sliceTriq = round.(UInt8, median_filter(sliceTriq)) + end + sliceTriq = reverse(sliceTriq, dims=2) + save_bitmap(joinpath("public", "TrIQ_$(text_nmass).bmp"), sliceTriq, ViridisPalette) + + imgIntT = "/TrIQ_$(text_nmass).bmp?t=$(timestamp)" + plotdataImgT, plotlayoutImgT, imgWidth, imgHeight = loadImgPlot(imgIntT) + current_triq = "TrIQ_$(text_nmass).bmp" + msgtriq = "TrIQ image with the Nmass of $(replace(text_nmass, "_" => "."))" + + bound = MSI_src.get_outlier_thres(slice, triqProb) + levels = range(bound[1], stop=bound[2], length=8) + levels = vcat(levels, 2 * levels[end] - levels[end-1]) + Colorbar(fig[1, 1], colormap=cgrad(:viridis, colorLevel, categorical=true), limits=(0, bound[2]), ticks=levels, tickformat=log_tick_formatter, label="Intensity", size=25) + save("public/colorbar_TrIQ_$(text_nmass).png", fig) + colorbarT = "/colorbar_TrIQ_$(text_nmass).png?t=$(timestamp)" + current_col_triq = "colorbar_TrIQ_$(text_nmass).png" + + triq_bmp = sort(filter(filename -> startswith(filename, "TrIQ_") && endswith(filename, ".bmp"), readdir("public")), lt=natural) + col_triq_png = sort(filter(filename -> startswith(filename, "colorbar_TrIQ_") && endswith(filename, ".png"), readdir("public")), lt=natural) + + fTime = time() + eTime = round(fTime - sTime, digits=3) + msg = "The TrIQ image has been created in $(eTime) seconds successfully inside the 'public' folder of the app" + selectedTab = "tab1" + end + else # If we don't use TrIQ + sliceQuant = quantize_intensity(slice, colorLevel) + if MFilterEnabled + sliceQuant = round.(UInt8, median_filter(sliceQuant)) + end + sliceQuant = reverse(sliceQuant, dims=2) + save_bitmap(joinpath("public", "MSI_$(text_nmass).bmp"), sliceQuant, ViridisPalette) + + imgInt = "/MSI_$(text_nmass).bmp?t=$(timestamp)" + plotdataImg, plotlayoutImg, imgWidth, imgHeight = loadImgPlot(imgInt) + current_msi = "MSI_$(text_nmass).bmp" + msgimg = "Image with the Nmass of $(replace(text_nmass, "_" => "."))" + + levels = range(0, maximum(slice), length=8) + Colorbar(fig[1, 1], colormap=cgrad(:viridis, colorLevel, categorical=true), limits=(0, maximum(slice)), ticks=levels, tickformat=log_tick_formatter, label="Intensity", size=25) + save("public/colorbar_MSI_$(text_nmass).png", fig) + colorbar = "/colorbar_MSI_$(text_nmass).png?t=$(timestamp)" + current_col_msi = "colorbar_MSI_$(text_nmass).png" + + msi_bmp = sort(filter(filename -> startswith(filename, "MSI_") && endswith(filename, ".bmp"), readdir("public")), lt=natural) + col_msi_png = sort(filter(filename -> startswith(filename, "colorbar_MSI_") && endswith(filename, ".png"), readdir("public")), lt=natural) + selectedTab = "tab0" + fTime = time() + eTime = round(fTime - sTime, digits=3) + msg = "The image has been created in $(eTime) seconds successfully inside the 'public' folder of the app" end catch e - msg="Warning: $e" + msg = "There was an error creating the image: $e" + warning_msg = true + @error "Image creation failed" exception=(e, catch_backtrace()) end - sliceQuant=reverse(sliceQuant, dims=2) - SaveBitmapCl(joinpath("public", "MSI_$(text_nmass).bmp"),sliceQuant,ViridisPalette) - # Use timestamp to refresh image interface container - imgInt="/MSI_$(text_nmass).bmp?t=$(timestamp)" - plotdataImg, plotlayoutImg, imgWidth, imgHeight=loadImgPlot(imgInt) - # Get current image - current_msi="MSI_$(text_nmass).bmp" - msgimg="Image with the Nmass of $(replace(text_nmass, "_" => "."))" - # Create colorbar - levels=range(0,maximum(slice),length=8) - Colorbar(fig[1, 1], colormap=cgrad(:viridis, colorLevel, categorical=true), limits=(0, maximum(slice)),ticks=levels,tickformat=log_tick_formatter, label="Intensity", size=25) - save("public/colorbar_MSI_$(text_nmass).png", fig) - colorbar="/colorbar_MSI_$(text_nmass).png?t=$(timestamp)" - # Get current colorbar - current_col_msi="colorbar_MSI_$(text_nmass).png" - # We update the directory to include the new placed images. - msi_bmp=sort(filter(filename -> startswith(filename, "MSI_") && endswith(filename, ".bmp"), readdir("public")),lt=natural) - col_msi_png=sort(filter(filename -> startswith(filename, "colorbar_MSI_") && endswith(filename, ".png"), readdir("public")),lt=natural) - selectedTab="tab0" - fTime=time() - eTime=round(fTime-sTime,digits=3) - msg="The file has been created in $(eTime) seconds successfully inside the 'public' folder of the app" + else + msg = "Invalid parameters. Nmass, Tol, or colorLevel are incorrect." + warning_msg = true + @error msg end - catch e - msg="There was an error loading the ImzML file, please verify the file accordingly and try again. $(e)" - warning_msg=true + finally + # This block will always run at the end of the async task + GC.gc() + if Sys.islinux() + ccall(:malloc_trim, Int32, (Int32,), 0) + end + btnStartDisable = false + btnPlotDisable = false + btnOpticalDisable = false + progress = false + btnSpectraDisable = false + SpectraEnabled = true end - else - msg="File does not exist or a parameter is incorrect, please try again." - warning_msg=true - end - GC.gc() # Trigger garbage collection - if Sys.islinux() - ccall(:malloc_trim, Int32, (Int32,), 0) # Ensure julia returns the freed memory to OS - end - btnStartDisable=false - btnPlotDisable=false - btnOpticalDisable=false - progress=false - if isfile(full_routeMz) - # We enable coord search and spectra plot creation - btnSpectraDisable=false - SpectraEnabled=true end end @onbutton createMeanPlot begin - msg="Mean spectrum plot selected" - sTime=time() - if isfile(full_routeMz) # Check if the file exists - progressSpectraPlot=true - btnPlotDisable=true - btnStartDisable=true - msg="Loading plot..." - plotdata, plotlayout, xSpectraMz, ySpectraMz=meanSpectrumPlot(full_routeMz) - GC.gc() # Trigger garbage collection - if Sys.islinux() - ccall(:malloc_trim, Int32, (Int32,), 0) # Ensure julia returns the freed memory to OS + if msi_data === nothing + msg = "No data loaded. Please select a file first." + warning_msg = true + return + end + + # UI updates immediately + progressSpectraPlot = true + btnPlotDisable = true + btnStartDisable = true + msg = "Loading plot..." + + @async begin + try + sTime = time() + plotdata, plotlayout, xSpectraMz, ySpectraMz = meanSpectrumPlot(msi_data) + + selectedTab = "tab2" + fTime = time() + eTime = round(fTime - sTime, digits=3) + msg = "Plot loaded in $(eTime) seconds" + catch e + msg = "Could not generate mean spectrum plot: $e" + warning_msg = true + @error "Mean spectrum plotting failed" exception=(e, catch_backtrace()) + finally + # This runs after the async task is finished + progressSpectraPlot = false + btnPlotDisable = false + btnSpectraDisable = false + if msi_data !== nothing && msi_data.source isa ImzMLSource + btnStartDisable = false + end end - selectedTab="tab2" - fTime=time() - eTime=round(fTime-sTime,digits=3) - msg="Plot loaded in $(eTime) seconds" - else - msg="there was an error with the mzML, please try again" - warning_msg=true end - progressSpectraPlot=false - btnPlotDisable=false - if endswith(full_route, "imzML") - btnStartDisable=false + end + + @onbutton createSumPlot begin + if msi_data === nothing + msg = "No data loaded. Please select a file first." + warning_msg = true + return end - if isfile(full_routeMz) - # We enable coord search and spectra plot creation - btnSpectraDisable=false - SpectraEnabled=true + + # UI updates immediately + progressSpectraPlot = true + btnPlotDisable = true + btnStartDisable = true + msg = "Loading total spectrum plot..." + + @async begin + try + sTime = time() + plotdata, plotlayout, xSpectraMz, ySpectraMz = sumSpectrumPlot(msi_data) + + selectedTab = "tab2" + fTime = time() + eTime = round(fTime - sTime, digits=3) + msg = "Total plot loaded in $(eTime) seconds" + catch e + msg = "Could not generate total spectrum plot: $e" + warning_msg = true + @error "Total spectrum plotting failed" exception=(e, catch_backtrace()) + finally + # This runs after the async task is finished + progressSpectraPlot = false + btnPlotDisable = false + btnSpectraDisable = false + if msi_data !== nothing && msi_data.source isa ImzMLSource + btnStartDisable = false + end + end end end @onbutton createXYPlot begin - msg="XY spectrum plot selected" - sTime=time() - if isfile(full_routeMz) # Check if the file exists - progressSpectraPlot=true - btnStartDisable=true - btnPlotDisable=true - btnSpectraDisable=true - msg="Loading plot..." - spectraMz=LoadMzml(full_routeMz) - layoutSpectra=PlotlyBase.Layout( - title="($xCoord, $yCoord) Specific spectrum plot", - xaxis=PlotlyBase.attr( - title="m/z", - showgrid=true - ), - yaxis=PlotlyBase.attr( - title="Intensity", - showgrid=true - ), - autosize=false, - margin=attr(l=0,r=0,t=120,b=0,pad=0) - ) - if xCoord < 1 - xCoord=1 - elseif xCoord > imgWidth - xCoord=imgWidth - end - if yCoord > -1 - yCoord=-1 - elseif yCoord < -imgHeight - yCoord=-imgHeight - end - xSpectraMz=spectraMz[1,abs(xCoord)] - ySpectraMz=spectraMz[2,abs(yCoord)] - traceSpectra=PlotlyBase.scatter(x=xSpectraMz, y=ySpectraMz, mode="lines") - plotdata=[traceSpectra] # We add the data from spectra to the plot - plotlayout=layoutSpectra - GC.gc() # Trigger garbage collection - if Sys.islinux() - ccall(:malloc_trim, Int32, (Int32,), 0) # Ensure julia returns the freed memory to OS - end - selectedTab="tab2" - fTime=time() - eTime=round(fTime-sTime,digits=3) - msg="Plot loaded in $(eTime) seconds" - else - msg="there was an error with the mzML or the coordenates, please try again" - warning_msg=true + if msi_data === nothing + msg = "No data loaded. Please select a file first." + warning_msg = true + return end - progressSpectraPlot=false - btnPlotDisable=false - if endswith(full_route, "imzML") - btnStartDisable=false - end - if isfile(full_routeMz) - # We enable coord search and spectra plot creation - btnSpectraDisable=false - SpectraEnabled=true + + # UI updates immediately + progressSpectraPlot = true + btnStartDisable = true + btnPlotDisable = true + btnSpectraDisable = true + msg = "Loading plot..." + + @async begin + try + sTime = time() + # The UI uses negative Y values, so we adjust before calling the plot function + y = yCoord < 0 ? abs(yCoord) : yCoord + plotdata, plotlayout, xSpectraMz, ySpectraMz = xySpectrumPlot(msi_data, xCoord, y, imgWidth, imgHeight) + + # Update UI coordinates + xCoord = plotlayout.title == "Spectrum #$(xCoord)" ? xCoord : clamp(xCoord, 1, imgWidth) + yCoord = plotlayout.title == "Spectrum #$(xCoord)" ? 0 : -clamp(y, 1, imgHeight) + + selectedTab = "tab2" + + fTime = time() + eTime = round(fTime - sTime, digits=3) + msg = "Plot loaded in $(eTime) seconds" + catch e + msg = "Could not retrieve spectrum: $e" + warning_msg = true + @error "Spectrum plotting failed" exception=(e, catch_backtrace()) + finally + # This runs after the async task is finished + progressSpectraPlot = false + btnPlotDisable = false + btnSpectraDisable = false + if msi_data !== nothing && msi_data.source isa ImzMLSource + btnStartDisable = false + end + end end end @@ -796,14 +792,21 @@ include("./julia_imzML_visual.jl") cleaned_imgInt=replace(imgInt, r"\?.*" => "") cleaned_imgInt=lstrip(cleaned_imgInt, '/') var=joinpath( "./public", cleaned_imgInt ) - sTime=time() - if isfile(var) - progressPlot=true - btnPlotDisable=true - btnStartDisable=true - btnSpectraDisable=true + if !isfile(var) + msg="Image could not be 3d plotted" + warning_msg=true + return + end + + progressPlot=true + btnPlotDisable=true + btnStartDisable=true + btnSpectraDisable=true + + @async begin try + sTime=time() plotdata3d, plotlayout3d=loadSurfacePlot(imgInt) GC.gc() # Trigger garbage collection if Sys.islinux() @@ -816,34 +819,39 @@ include("./julia_imzML_visual.jl") catch e msg="Failed to load and process image: $e" warning_msg=true + finally + progressPlot=false + btnPlotDisable=false + btnStartDisable=false + if msi_data !== nothing + # We enable coord search and spectra plot creation + btnSpectraDisable=false + SpectraEnabled=true + end end - else - msg="Image could not be 3d plotted" - warning_msg=true end - progressPlot=false - btnPlotDisable=false - btnStartDisable=false - if isfile(full_routeMz) - # We enable coord search and spectra plot creation - btnSpectraDisable=false - SpectraEnabled=true - end - end - # 3d plot for TrIQ + end # 3d plot for TrIQ + @onbutton triq3dPlot begin msg="TrIQ 3D plot selected" cleaned_imgIntT=replace(imgIntT, r"\?.*" => "") cleaned_imgIntT=lstrip(cleaned_imgIntT, '/') var=joinpath( "./public", cleaned_imgIntT ) - sTime=time() - if isfile(var) - progressPlot=true - btnPlotDisable=true - btnStartDisable=true - btnSpectraDisable=true + if !isfile(var) + msg="Image could not be 3d plotted" + warning_msg=true + return + end + + progressPlot=true + btnPlotDisable=true + btnStartDisable=true + btnSpectraDisable=true + + @async begin try + sTime=time() plotdata3d, plotlayout3d=loadSurfacePlot(imgIntT) GC.gc() # Trigger garbage collection if Sys.islinux() @@ -856,18 +864,16 @@ include("./julia_imzML_visual.jl") catch e msg="Failed to load and process image: $e" warning_msg=true + finally + progressPlot=false + btnPlotDisable=false + btnStartDisable=false + if msi_data !== nothing + # We enable coord search and spectra plot creation + btnSpectraDisable=false + SpectraEnabled=true + end end - else - msg="Image could not be 3d plotted" - warning_msg=true - end - progressPlot=false - btnPlotDisable=false - btnStartDisable=false - if isfile(full_routeMz) - # We enable coord search and spectra plot creation - btnSpectraDisable=false - SpectraEnabled=true end end @@ -877,15 +883,21 @@ include("./julia_imzML_visual.jl") cleaned_imgInt=replace(imgInt, r"\?.*" => "") cleaned_imgInt=lstrip(cleaned_imgInt, '/') var=joinpath("./public", cleaned_imgInt) - sTime=time() - - if isfile(var) - progressPlot=true - btnPlotDisable=true - btnStartDisable=true - btnSpectraDisable=true + + if !isfile(var) + msg="Image could not be 2D plotted" + warning_msg=true + return + end + + progressPlot=true + btnPlotDisable=true + btnStartDisable=true + btnSpectraDisable=true + + @async begin try - img=load(var) + sTime=time() plotdataC,plotlayoutC=loadContourPlot(imgInt) GC.gc() # Trigger garbage collection if Sys.islinux() @@ -898,18 +910,16 @@ include("./julia_imzML_visual.jl") catch e msg="Failed to load and process image: $e" warning_msg=true + finally + progressPlot=false + btnPlotDisable=false + btnStartDisable=false + if msi_data !== nothing + # We enable coord search and spectra plot creation + btnSpectraDisable=false + SpectraEnabled=true + end end - else - msg="Image could not be 2D plotted" - warning_msg=true - end - progressPlot=false - btnPlotDisable=false - btnStartDisable=false - if isfile(full_routeMz) - # We enable coord search and spectra plot creation - btnSpectraDisable=false - SpectraEnabled=true end end # Contour 2d plot for TrIQ @@ -918,15 +928,21 @@ include("./julia_imzML_visual.jl") cleaned_imgIntT=replace(imgIntT, r"\?.*" => "") cleaned_imgIntT=lstrip(cleaned_imgIntT, '/') var=joinpath("./public", cleaned_imgIntT) - sTime=time() - - if isfile(var) - progressPlot=true - btnPlotDisable=true - btnStartDisable=true - btnSpectraDisable=true + + if !isfile(var) + msg="Image could not be 2D plotted" + warning_msg=true + return + end + + progressPlot=true + btnPlotDisable=true + btnStartDisable=true + btnSpectraDisable=true + + @async begin try - img=load(var) + sTime=time() plotdataC,plotlayoutC=loadContourPlot(imgIntT) GC.gc() # Trigger garbage collection if Sys.islinux() @@ -939,18 +955,16 @@ include("./julia_imzML_visual.jl") catch e msg="Failed to load and process image: $e" warning_msg=true + finally + progressPlot=false + btnPlotDisable=false + btnStartDisable=false + if msi_data !== nothing + # We enable coord search and spectra plot creation + btnSpectraDisable=false + SpectraEnabled=true + end end - else - msg="Image could not be 2D plotted" - warning_msg=true - end - progressPlot=false - btnPlotDisable=false - btnStartDisable=false - if isfile(full_routeMz) - # We enable coord search and spectra plot creation - btnSpectraDisable=false - SpectraEnabled=true end end @@ -961,48 +975,39 @@ include("./julia_imzML_visual.jl") # To include a visualization in the spectrum plot indicating where is the selected mass @onchange Nmass begin if !isempty(xSpectraMz) - traceSpectra=PlotlyBase.scatter(x=xSpectraMz, y=ySpectraMz, mode="lines",name="Spectra",showlegend=false) - trace2=PlotlyBase.scatter(x=[Nmass, Nmass],y=[0, maximum(ySpectraMz)],mode="lines",line=attr(color="red", width=0.5),name="m/z selected",showlegend=false) - plotdata=[traceSpectra,trace2] # We add the data from spectra and the red line to the plot + # Use a stem plot for the main spectrum for consistency + traceSpectra = PlotlyBase.stem(x=xSpectraMz, y=ySpectraMz, marker=attr(size=1, color="blue", opacity=0.5), name="Spectrum", hoverinfo="x", hovertemplate="m/z: %{x:.4f}", showlegend=false) + + # Keep this as a scatter plot to draw the vertical line + trace2 = PlotlyBase.scatter(x=[Nmass, Nmass], y=[0, maximum(ySpectraMz)], mode="lines", line=attr(color="red", width=0.5), name="m/z selected", showlegend=false) + + plotdata = [traceSpectra, trace2] end end # Event detection for clicking on the images @onchange data_click begin + if selectedTab == "tab1" || selectedTab == "tab0" + # This is for the image heatmaps + cursor_data = data_click["cursor"] + x = Int32(round(cursor_data["x"])) + y = Int32(round(cursor_data["y"])) # y is negative in the UI + + # Update the reactive coordinates, which will trigger the crosshair update + xCoord = clamp(x, 1, imgWidth) + yCoord = clamp(y, -imgHeight, -1) + end + end + + @onchange xCoord, yCoord begin if selectedTab == "tab1" - cursor_data=data_click["cursor"] - xCoord=Int32(round(cursor_data["x"])) - yCoord=Int32(round(cursor_data["y"])) - if xCoord < 1 - xCoord=1 - elseif xCoord > imgWidth - xCoord=imgWidth - end - if yCoord > -1 - yCoord=-1 - elseif yCoord < -imgHeight - yCoord=-imgHeight - end # Get the x and y values from the click of the cursor and make sure they don't exceed image proportions - plotdataImgT=filter(trace -> !(get(trace, :name, "") in ["Line X", "Line Y"]), plotdataImgT) - trace1, trace2=crossLinesPlot(xCoord, yCoord, imgWidth, -imgHeight) - plotdataImgT=append!(plotdataImgT, [trace1, trace2]) + plotdataImgT = filter(trace -> !(get(trace, :name, "") in ["Line X", "Line Y"]), plotdataImgT) + trace1, trace2 = crossLinesPlot(xCoord, yCoord, imgWidth, -imgHeight) + plotdataImgT = append!(plotdataImgT, [trace1, trace2]) elseif selectedTab == "tab0" - cursor_data=data_click["cursor"] - xCoord=Int32(round(cursor_data["x"])) - yCoord=Int32(round(cursor_data["y"])) - if xCoord < 1 - xCoord=1 - elseif xCoord > imgWidth - xCoord=imgWidth - end - if yCoord > -1 - yCoord=-1 - elseif yCoord < -imgHeight - yCoord=-imgHeight - end # Get the x and y values from the click of the cursor and make sure they don't exceed image proportions - plotdataImg=filter(trace -> !(get(trace, :name, "") in ["Line X", "Line Y","Optical"]), plotdataImg) - trace1, trace2=crossLinesPlot(xCoord, yCoord, imgWidth, -imgHeight) - plotdataImg=append!(plotdataImg, [trace1, trace2]) + plotdataImg = filter(trace -> !(get(trace, :name, "") in ["Line X", "Line Y", "Optical"]), plotdataImg) + trace1, trace2 = crossLinesPlot(xCoord, yCoord, imgWidth, -imgHeight) + plotdataImg = append!(plotdataImg, [trace1, trace2]) end end diff --git a/app.jl.html b/app.jl.html index 77c2180..f9c9653 100644 --- a/app.jl.html +++ b/app.jl.html @@ -69,6 +69,11 @@ Mean spectrum plot + + + Sum Spectrum plot + + Spectrum plot (X,Y) diff --git a/julia_imzML_visual.jl b/julia_imzML_visual.jl index 604fcb1..e20abf7 100644 --- a/julia_imzML_visual.jl +++ b/julia_imzML_visual.jl @@ -1,84 +1,3 @@ -# IntQuantCl is originally a function of julia mzMl imzML with the adition -# of altering the scale of the colors according to colorlevel. -function IntQuantCl( slice , colorLevel) - # Compute scale factor for amplitude discretization - lower = minimum( slice ) - scale = colorLevel / maximum( slice ) - dim = size( slice ) - image = zeros( UInt8, dim[1], dim[2] ) - for i in 1:length( slice ) - image[i] = convert( UInt8, floor( slice[i] * scale + 0.5 ) ) - end - return image -end - -# SaveBitmap originally a function of the mzML imzML library in julia, -# This function dinamically adjust the color palete adjusting to the ammount of colors -# available in pixmap -function SaveBitmapCl( name, pixMap::Array{UInt8,2}, colorTable::Array{UInt32,1} ) - # Get image dimensions - dim = size( pixMap ) - if length( dim ) != 2 - return 0 - end - # Normalize pixel values to get a more accurate reading of the image - minVal = minimum(pixMap) - maxVal = maximum(pixMap) - pixMap = round.(UInt8, 255 * (pixMap .- minVal) ./ (maxVal - minVal)) - # Compute row padding - padding = ( 4 - dim[1] & 0x3 ) & 0x3 - # Compute file dimensions. Header = 14 + 40 + ( 256 * 4 ) = 1078 - offset = 1078 - imgBytes = dim[2] * ( dim[1] + padding ) - # Create file - stream = open( name, "w" ) - # Save file header - write( stream, UInt16( 0x4D42 ) ) - write( stream, UInt32[ offset + imgBytes, 0 , offset ] ) - # Save info header - write( stream, UInt32[ 40, dim[1], dim[2], 0x80001, 0 ] ) - write( stream, UInt32[ imgBytes, 0, 0, 256, 0 ] ) - # Save color table - write( stream, colorTable ) - if length( colorTable ) < 256 - fixTable = zeros( UInt32, 256 - length( colorTable ) ) - write( stream, fixTable ) - end - # Save image pixels - if padding == 0 - for i = 1:dim[2] - write( stream, pixMap[:,i] ) - end - else - zeroPad = zeros( UInt8, padding ) - for i in 1:dim[2] - write( stream, pixMap[:,i] ) - write( stream, zeroPad ) - end - end - # Close file - close( stream ) -end - -# SaveBitmap originally a function of the mzML imzML library in julia, -# now has an adjustment for NaN values in case they exist to mantain data integrity -function GetMzSliceJl(imzML, mass, tolerance) - # Alloc space for slice - width = maximum(imzML[1, :]) - height = maximum(imzML[2, :]) - image = fill(0.0, width, height) - - for i in 1:size(imzML)[2] - index = julia_mzML_imzML.FindMass(imzML[3, i], mass, tolerance) - if index != 0 - image[imzML[1, i], imzML[2, i]] = imzML[4, i][index] - end - end - # Adjustment for NaN values with 0 - replace!(image, NaN => 0.0) - return image -end - # == Search functions == # Functions that recieve a list to update, and the current direction both as string for # searching in the directory the position the list is going @@ -129,10 +48,12 @@ function loadImgPlot(interfaceImg::String) layout=PlotlyBase.Layout( xaxis=PlotlyBase.attr( visible=false, - scaleanchor="y" + scaleanchor="y", + range=[0, width] ), yaxis=PlotlyBase.attr( - visible=false + visible=false, + range=[-height, 0] ), margin=attr(l=0,r=0,t=0,b=0,pad=0) ) @@ -167,6 +88,7 @@ function loadImgPlot(interfaceImg::String) plotlayout=layout return plotdata, plotlayout, width, height end + # loadImgPlot recieves the local directory of the image as a string, the local directory o the overlay image # and the transparency its required to have. Returns the layout and data for the heatmap plotly plot # this function loads the image into a plot @@ -286,6 +208,7 @@ function loadContourPlot(interfaceImg::String) plotlayout=layout return plotdata, plotlayout end + # loadSurfacePlot recieves the local directory of the image as a string, # returns the layout and data for the surface plotly plot # this function loads the image and applies a gaussian filter @@ -363,6 +286,7 @@ function loadSurfacePlot(interfaceImg::String) plotlayout=layout3D return plotdata, plotlayout end + # This function recieves the x and y coords currently selected, and the dimentions of # the image to create two traces that will display in a cross section function crossLinesPlot(x, y, maxwidth, maxheight) @@ -404,77 +328,117 @@ function log_tick_formatter(values::Vector{Float64}) end -# Median filter: an adaptation of the R medianfilter, which averages the matrix -# with the close pixels just from the sides to reduce noise. -# This one in particular is a midpoint fiter from a 3x3 neighbour area -function medianFilterjl(pixMap) - height, width = size(pixMap) - padded_pixMap = padarray(pixMap, 1) - target = zeros(eltype(pixMap), height, width) - - for j in 1:width - for i in 1:height - neighbors = [] - for dj in j:j+2 - for di in i:i+2 - push!(neighbors, padded_pixMap[di, dj]) - end - end - target[i, j] = median(neighbors) - end - end - - return target -end - -function padarray(A, padsize) - h, w = size(A) - padded = zeros(eltype(A), h + 2*padsize, w + 2*padsize) - padded[padsize+1:end-padsize, padsize+1:end-padsize] .= A - padded[1:padsize, padsize+1:end-padsize] .= A[1:padsize, :] - padded[end-padsize+1:end, padsize+1:end-padsize] .= A[end-padsize+1:end, :] - padded[padsize+1:end-padsize, 1:padsize] .= A[:, 1:padsize] - padded[padsize+1:end-padsize, end-padsize+1:end] .= A[:, end-padsize+1:end] - padded[1:padsize, 1:padsize] .= A[1, 1] - padded[1:padsize, end-padsize+1:end] .= A[1, end] - padded[end-padsize+1:end, 1:padsize] .= A[end, 1] - padded[end-padsize+1:end, end-padsize+1:end] .= A[end, end] - return padded -end - # meanSpectrumPlot recieves the local directory of the image as a string, # returns the layout and data for the surface plotly plot # this function loads the spectra data and makes a mean to display # its values in the spectrum plot -function meanSpectrumPlot(mzmlRoute::String) - spectraMz=LoadMzml(mzmlRoute) - xSpectraMz=Float64[] - ySpectraMz=Float64[] - - layout=PlotlyBase.Layout( - title="Mean spectrum plot", +function meanSpectrumPlot(data::MSIData) + layout = PlotlyBase.Layout( + title="Average Spectrum Plot", + hovermode="closest", xaxis=PlotlyBase.attr( title="m/z", + showgrid=true + ), + yaxis=PlotlyBase.attr( + title="Average Intensity", showgrid=true, - tickformat = ".3g" + tickformat=".3g" + ), + margin=attr(l=0, r=0, t=120, b=0, pad=0) + ) + + # Use the new, efficient function from the backend + xSpectraMz, ySpectraMz = get_average_spectrum(data) + + if isempty(xSpectraMz) + @warn "Average spectrum is empty." + trace = PlotlyBase.stem(x=Float64[], y=Float64[]) + else + trace = PlotlyBase.stem(x=xSpectraMz, y=ySpectraMz, marker=attr(size=1, color="blue", opacity=0.5), name="Average", hoverinfo="x",hovertemplate="m/z: %{x:.4f}") + end + + + plotdata = [trace] + plotlayout = layout + return plotdata, plotlayout, xSpectraMz, ySpectraMz +end + +function xySpectrumPlot(data::MSIData, xCoord::Int, yCoord::Int, imgWidth::Int, imgHeight::Int) + local mz::AbstractVector, intensity::AbstractVector + local plot_title::String + + is_imaging = data.source isa ImzMLSource + + if is_imaging + # For imaging data, use (X, Y) coordinates + x = clamp(xCoord, 1, imgWidth) + y = clamp(yCoord, 1, imgHeight) + + mz, intensity = GetSpectrum(data, x, y) + plot_title = "Spectrum at ($x, $y)" + else + # For non-imaging data, treat xCoord as the spectrum index + index = clamp(xCoord, 1, length(data.spectra_metadata)) + + mz, intensity = GetSpectrum(data, index) + plot_title = "Spectrum #$index" + end + + layout = PlotlyBase.Layout( + title=plot_title, + hovermode="closest", + xaxis=PlotlyBase.attr( + title="m/z", + showgrid=true ), yaxis=PlotlyBase.attr( title="Intensity", showgrid=true, - tickformat = ".3g" + tickformat=".3g" ), - autosize=false, - margin=attr(l=0,r=0,t=120,b=0,pad=0) + margin=attr(l=0, r=0, t=120, b=0, pad=0) ) - try - xSpectraMz=mean(spectraMz[1,:]) - ySpectraMz=mean(spectraMz[2,:]) - catch e - xSpectraMz=spectraMz[1,1] - ySpectraMz=spectraMz[2,1] + + # Downsample for plotting performance + mz_down, int_down = MSI_src.downsample_spectrum(mz, intensity) + + trace = PlotlyBase.stem(x=mz_down, y=int_down, marker=attr(size=1, color="blue", opacity=0.5), name="Spectrum", hoverinfo="x", hovertemplate="m/z: %{x:.4f}") + + plotdata = [trace] + plotlayout = layout + + # Return the full data for other uses, and the plot data + return plotdata, plotlayout, mz, intensity +end + +function sumSpectrumPlot(data::MSIData) + layout = PlotlyBase.Layout( + title="Total Spectrum Plot", + hovermode="closest", + xaxis=PlotlyBase.attr( + title="m/z", + showgrid=true + ), + yaxis=PlotlyBase.attr( + title="Total Intensity", + showgrid=true, + tickformat=".3g" + ), + margin=attr(l=0, r=0, t=120, b=0, pad=0) + ) + + # Use the get_total_spectrum function from the backend + xSpectraMz, ySpectraMz = get_total_spectrum(data) + + if isempty(xSpectraMz) + @warn "Total spectrum is empty." + trace = PlotlyBase.stem(x=Float64[], y=Float64[]) + else + trace = PlotlyBase.stem(x=xSpectraMz, y=ySpectraMz, marker=attr(size=1, color="blue", opacity=0.5), name="Total", hoverinfo="x",hovertemplate="m/z: %{x:.4f}") end - trace=PlotlyBase.stem(x=xSpectraMz, y=ySpectraMz,marker=attr(size=1, color="blue", opacity=0.1)) - plotdata=[trace] # We add the data from spectra to the plot - plotlayout=layout + + plotdata = [trace] + plotlayout = layout return plotdata, plotlayout, xSpectraMz, ySpectraMz -end \ No newline at end of file +end diff --git a/src/MSIData.jl b/src/MSIData.jl index 4ac7107..b570ee3 100644 --- a/src/MSIData.jl +++ b/src/MSIData.jl @@ -6,7 +6,7 @@ including caching and iteration logic, for handling large mzML and imzML dataset efficiently. """ -using Base64, Libz # For reading binary data +using Base64, Libz, Serialization # For reading binary data # Abstract type for different data sources (e.g., mzML, imzML) # This allows dispatching to the correct binary reading logic. @@ -179,17 +179,10 @@ function GetSpectrum(data::MSIData, x::Int, y::Int) return GetSpectrum(data, index) # Call the existing method end -""" - get_total_spectrum(msi_data::MSIData; num_bins::Int=20000) -> Tuple{Vector{Float64}, Vector{Float64}} +using Serialization -Calculates the sum of all spectra in the dataset by binning. -The m/z range is determined dynamically by finding the global min/max m/z values -across the entire dataset. - -Returns a tuple containing two vectors: the binned m/z axis and the summed intensities. -""" -function get_total_spectrum(msi_data::MSIData; num_bins::Int=20000) - println("Calculating total spectrum...") +function get_total_spectrum_imzml(msi_data::MSIData; num_bins::Int=2000) + println("Calculating total spectrum (2-pass method for imzML)...") # 1. First Pass: Find the global m/z range println(" Pass 1: Finding global m/z range...") @@ -222,7 +215,6 @@ function get_total_spectrum(msi_data::MSIData; num_bins::Int=20000) return # continue equivalent end for i in eachindex(mz) - # Find the correct bin for the current m/z value bin_index = clamp(round(Int, (mz[i] - global_min_mz) / bin_step) + 1, 1, num_bins) intensity_sum[bin_index] += intensity[i] end @@ -232,6 +224,83 @@ function get_total_spectrum(msi_data::MSIData; num_bins::Int=20000) return (collect(mz_bins), intensity_sum) end +function get_total_spectrum_mzml(msi_data::MSIData; num_bins::Int=2000) + println("Calculating total spectrum (single-pass optimization for mzML)...") + num_spectra = length(msi_data.spectra_metadata) + if num_spectra == 0 + return (Float64[], Float64[]) + end + + temp_path, temp_io = mktemp() + try + # --- Pass 1: Write spectra to temp file and find min/max m/z --- + println(" Pass 1: Caching spectra and finding global m/z range...") + global_min_mz = Inf + global_max_mz = -Inf + + _iterate_spectra_fast(msi_data) do idx, mz, intensity + Serialization.serialize(temp_io, (mz, intensity)) + if !isempty(mz) + local_min, local_max = extrema(mz) + global_min_mz = min(global_min_mz, local_min) + global_max_mz = max(global_max_mz, local_max) + end + end + + flush(temp_io) + + if !isfinite(global_min_mz) + @warn "Could not determine a valid m/z range. All spectra might be empty." + return (Float64[], Float64[]) + end + println(" Global m/z range found: [$(global_min_mz), $(global_max_mz)]") + + # --- Pass 2: Read from temp file and bin intensities --- + println(" Pass 2: Reading from cache and summing intensities into $num_bins bins...") + + seekstart(temp_io) + mz_bins = range(global_min_mz, stop=global_max_mz, length=num_bins) + intensity_sum = zeros(Float64, num_bins) + bin_step = step(mz_bins) + + while !eof(temp_io) + mz, intensity = Serialization.deserialize(temp_io)::Tuple{AbstractVector, AbstractVector} + + if isempty(mz) + continue + end + for i in eachindex(mz) + bin_index = clamp(round(Int, (mz[i] - global_min_mz) / bin_step) + 1, 1, num_bins) + intensity_sum[bin_index] += intensity[i] + end + end + + println("Total spectrum calculation complete.") + return (collect(mz_bins), intensity_sum) + + finally + close(temp_io) + rm(temp_path, force=true) + end +end + +""" + get_total_spectrum(msi_data::MSIData; num_bins::Int=20000) -> Tuple{Vector{Float64}, Vector{Float64}} + +Calculates the sum of all spectra in the dataset by binning. +This function dispatches to a specialized implementation based on the file type +(.imzML or .mzML) for optimal performance. + +Returns a tuple containing two vectors: the binned m/z axis and the summed intensities. +""" +function get_total_spectrum(msi_data::MSIData; num_bins::Int=2000) + if msi_data.source isa ImzMLSource + return get_total_spectrum_imzml(msi_data, num_bins=num_bins) + else # MzMLSource + return get_total_spectrum_mzml(msi_data, num_bins=num_bins) + end +end + """ get_average_spectrum(msi_data::MSIData; num_bins::Int=20000) -> Tuple{Vector{Float64}, Vector{Float64}} @@ -241,7 +310,7 @@ This is effectively the total ion chromatogram (TIC) divided by the number of sp Returns a tuple containing two vectors: the binned m/z axis and the averaged intensities. """ -function get_average_spectrum(msi_data::MSIData; num_bins::Int=20000) +function get_average_spectrum(msi_data::MSIData; num_bins::Int=2000) # This function uses the exact same logic as get_total_spectrum... mz_bins, intensity_sum = get_total_spectrum(msi_data, num_bins=num_bins) @@ -365,5 +434,3 @@ function _iterate_spectra_fast(f::Function, data::MSIData) # Dispatch to the correct implementation based on the source type _iterate_spectra_fast_impl(f, data, data.source) end - - diff --git a/src/MzmlConverter.jl b/src/MzmlConverter.jl index 69e10be..695eed4 100644 --- a/src/MzmlConverter.jl +++ b/src/MzmlConverter.jl @@ -26,26 +26,19 @@ struct BinaryMetadata end """ - GetMzmlScanTime(fileName::String) + GetMzmlScanTime_linebyline(fileName::String) -Parses a .mzML file to extract the scan start time for each spectrum. - -# Arguments -* `fileName`: Path to the .mzML file. - -# Returns -- A `Matrix{Int64}` where each row is `[spectrum_index, time_in_milliseconds]`. +(Internal) Slow, line-by-line parser for scan times. Used as a fallback if +the .mzML file index is missing or corrupt. """ -function GetMzmlScanTime(fileName::String) +function GetMzmlScanTime_linebyline(fileName::String) times = Tuple{Int64, Int64}[] - # Pre-allocate based on an estimate to reduce re-allocations try file_size = filesize(fileName) - estimated_spectra = max(1000, file_size ÷ 10000) # Heuristic + estimated_spectra = max(1000, file_size ÷ 10000) sizehint!(times, estimated_spectra) catch - # Ignore if filesize fails, sizehint! is just an optimization end open(fileName, "r") do stream @@ -60,27 +53,23 @@ function GetMzmlScanTime(fileName::String) current_index = 0 current_time = nothing - # Process index from the start tag itself idx_match = match(r"index=\"(\d+)\"", line) if idx_match !== nothing current_index = parse(Int, idx_match.captures[1]) + 1 end end elseif state == :in_spectrum - # Look for time cvParam if occursin("", line) state = :outside_spectrum if current_index > 0 && current_time !== nothing @@ -91,9 +80,100 @@ function GetMzmlScanTime(fileName::String) end end - # Sort by spectrum index to ensure correct order before normalization sort!(times, by = x -> x[1]) + if !isempty(times) + first_time = times[1][2] + for i in eachindex(times) + times[i] = (times[i][1], times[i][2] - first_time) + end + end + + if isempty(times) + return Matrix{Int64}(undef, 0, 2) + end + return permutedims(hcat(collect.(times)...)) +end + +""" + GetMzmlScanTime(fileName::String) + +Parses a .mzML file to extract the scan start time for each spectrum. +Uses the indexed part of the .mzML file for fast access, falling back to +a slower line-by-line parse if the index is not present. + +# Arguments +* `fileName`: Path to the .mzML file. + +# Returns +- A `Matrix{Int64}` where each row is `[spectrum_index, time_in_milliseconds]`. +""" +function GetMzmlScanTime(fileName::String) + times = Tuple{Int64, Int64}[] + + try + open(fileName, "r") do stream + # 1. Find and parse the spectrum index offsets + seekend(stream) + end_chunk_size = min(filesize(stream), 8192) + seek(stream, filesize(stream) - end_chunk_size) + footer = read(stream, String) + + index_offset_match = match(r"(\d+)", footer) + if index_offset_match === nothing + @warn "No found. Falling back to slow line-by-line parsing for scan times. This may be memory intensive." + return GetMzmlScanTime_linebyline(fileName) + end + + index_offset = parse(Int64, index_offset_match.captures[1]) + seek(stream, index_offset) + + # The find_tag function is defined in ParserHelpers.jl + if find_tag(stream, r"", line) || occursin("", line) + break + end + end + + if current_time !== nothing + push!(times, (idx, current_time)) + end + end + end + catch e + @error "Failed to parse scan times with indexed method. Falling back to line-by-line." exception=(e, catch_backtrace()) + return GetMzmlScanTime_linebyline(fileName) + end + + # The index is already sorted by spectrum index, so no need to sort `times`. + # Normalize times relative to the first scan if !isempty(times) first_time = times[1][2] @@ -102,7 +182,6 @@ function GetMzmlScanTime(fileName::String) end end - # Convert vector of tuples to a matrix if isempty(times) return Matrix{Int64}(undef, 0, 2) end @@ -329,131 +408,113 @@ function RenderPixel(pixel_info, scans, msi_data::MSIData, scan_time_deltas, pix return (mz_array, new_intensity) end -function ConvertMzmlToImzml(source_file::String, timing_matrix::Matrix{Int64}, scans::Matrix{Int64}) +function ConvertMzmlToImzml(source_file::String, target_ibd_file::String, timing_matrix::Matrix{Int64}, scans::Matrix{Int64}) if size(timing_matrix, 1) == 0 - return ProcessedPixel[], (0, 0) + # Create an empty .ibd file if there's nothing to process + open(target_ibd_file, "w") do ibd_stream + write(ibd_stream, zeros(UInt8, 16)) # UUID placeholder + end + return BinaryMetadata[], Tuple{Int, Int}[], (0, 0) end - # Get image dimensions from timing_matrix (assuming columns 1 and 2 are X and Y) width = maximum(timing_matrix[:, 1]) height = maximum(timing_matrix[:, 2]) - # Open the mzML file using the new memory-efficient API msi_data = OpenMSIData(source_file) - # Pre-calculate time deltas robustly scan_time_deltas = zeros(Int64, size(scans, 1)) if size(scans, 1) > 1 for i in 1:(size(scans, 1) - 1) delta = scans[i+1, 2] - scans[i, 2] - scan_time_deltas[i] = max(1, delta) # Ensure delta is at least 1 + scan_time_deltas[i] = max(1, delta) end - scan_time_deltas[end] = max(1, scan_time_deltas[end-1]) # Use previous delta for last scan, ensure at least 1 + scan_time_deltas[end] = max(1, scan_time_deltas[end-1]) end pixel_time_deltas = zeros(Int64, size(timing_matrix, 1)) if size(timing_matrix, 1) > 1 for i in 1:(size(timing_matrix, 1) - 1) - # Assumes time is in the 3rd column of the timing_matrix delta = timing_matrix[i+1, 3] - timing_matrix[i, 3] - pixel_time_deltas[i] = max(1, delta) # Ensure delta is at least 1 + pixel_time_deltas[i] = max(1, delta) end - pixel_time_deltas[end] = max(1, pixel_time_deltas[end-1]) # Use previous delta for last pixel, ensure at least 1 + pixel_time_deltas[end] = max(1, pixel_time_deltas[end-1]) end - processed_pixels = ProcessedPixel[] - sizehint!(processed_pixels, size(timing_matrix, 1)) + binary_meta_vec = BinaryMetadata[] + sizehint!(binary_meta_vec, size(timing_matrix, 1)) + coords_vec = Tuple{Int, Int}[] + sizehint!(coords_vec, size(timing_matrix, 1)) empty_pixel_count = 0 - # DEBUG: Scan coverage analysis - println("DEBUG: Scan coverage analysis:") - covered_pixels_count = 0 - for i in 1:size(timing_matrix, 1) - first_scan = timing_matrix[i, 4] - last_scan = timing_matrix[i, 5] - if first_scan <= last_scan && first_scan >= 1 && last_scan <= size(scans, 1) - covered_pixels_count += 1 + open(target_ibd_file, "w") do ibd_stream + write(ibd_stream, zeros(UInt8, 16)) # UUID placeholder + + for i in 1:size(timing_matrix, 1) + pixel_info = timing_matrix[i, :] + x, y = pixel_info[1], pixel_info[2] + push!(coords_vec, (x, y)) + + first_scan = pixel_info[4] + last_scan = pixel_info[5] + + if first_scan > last_scan || first_scan < 1 || last_scan > size(scans, 1) + empty_pixel_count += 1 + # For empty pixels, offsets point to the current end of file, with zero length + current_pos = position(ibd_stream) + push!(binary_meta_vec, BinaryMetadata(current_pos, 0, current_pos, 0)) + continue + end + + mz, intensity = RenderPixel(pixel_info, scans, msi_data, scan_time_deltas, pixel_time_deltas) + + # Write m/z array + mz_offset = position(ibd_stream) + for val in mz + write(ibd_stream, htol(Float64(val))) + end + mz_length = position(ibd_stream) - mz_offset + + # Write intensity array + int_offset = position(ibd_stream) + for val in intensity + write(ibd_stream, htol(Float32(val))) + end + int_length = position(ibd_stream) - int_offset + + push!(binary_meta_vec, BinaryMetadata(mz_offset, mz_length, int_offset, int_length)) end end - println(" Pixels with potential scans: $covered_pixels_count/$(size(timing_matrix, 1))") - println(" Total scans available: $(size(scans, 1))") - for i in 1:size(timing_matrix, 1) - pixel_info = timing_matrix[i, :] - x = pixel_info[1] - y = pixel_info[2] - first_scan = pixel_info[4] - last_scan = pixel_info[5] - - if first_scan > last_scan || first_scan < 1 || last_scan > size(scans, 1) - empty_pixel_count += 1 - push!(processed_pixels, ProcessedPixel((x, y), Float64[], Float32[])) - continue - end - - mz, intensity = RenderPixel(pixel_info, scans, msi_data, scan_time_deltas, pixel_time_deltas) - push!(processed_pixels, ProcessedPixel((x, y), mz, intensity)) - end - - @info "Found and created $empty_pixel_count empty pixels out of $(size(timing_matrix, 1)) total." - return processed_pixels, (width, height) + @info "Found and processed $empty_pixel_count empty pixels out of $(size(timing_matrix, 1)) total." + return binary_meta_vec, coords_vec, (width, height) end -function ExportImzml(target_file::String, pixels::Vector{ProcessedPixel}, dims::Tuple{Int, Int}) +function ExportImzml(target_file::String, binary_meta::Vector{BinaryMetadata}, coords::Vector{Tuple{Int, Int}}, dims::Tuple{Int, Int}) ibd_file = replace(target_file, r"\.imzML$"i => ".ibd") - # Ensure we have valid pixels - if isempty(pixels) - @error "No pixels to export" - return false + if isempty(binary_meta) + @warn "No binary metadata to export; creating empty imzML file." + # Still create a valid, empty imzML file end - binary_meta = BinaryMetadata[] - sizehint!(binary_meta, length(pixels)) - try - # Step 1: Write .ibd file with proper binary format - open(ibd_file, "w") do ibd_stream - # Write UUID as first 16 bytes (placeholder) - write(ibd_stream, zeros(UInt8, 16)) + # The .ibd file is now written by ConvertMzmlToImzml. + # This function is only responsible for the .imzML XML metadata file. - # Write in little-endian format (standard for imzML) - for pixel in pixels - # Write m/z array as 64-bit little-endian floats - mz_offset = position(ibd_stream) - for val in pixel.mz - write(ibd_stream, htol(Float64(val))) # Host to little-endian - end - mz_length = position(ibd_stream) - mz_offset - - # Write intensity array as 32-bit little-endian floats - int_offset = position(ibd_stream) - for val in pixel.intensity - write(ibd_stream, htol(Float32(val))) # Host to little-endian - end - int_length = position(ibd_stream) - int_offset - - push!(binary_meta, BinaryMetadata(mz_offset, mz_length, int_offset, int_length)) - end - end - - # Step 2: Write proper .imzML XML file open(target_file, "w") do imzml_stream - # Use indexedmzML wrapper + # XML Header write(imzml_stream, """ """) - # CV List + # CV List, File Description, etc. (static parts) write(imzml_stream, """ """) - - # File Description write(imzml_stream, """ @@ -462,8 +523,6 @@ function ExportImzml(target_file::String, pixels::Vector{ProcessedPixel}, dims:: """) - - # Referenceable Param Groups write(imzml_stream, """ @@ -479,24 +538,18 @@ function ExportImzml(target_file::String, pixels::Vector{ProcessedPixel}, dims:: """) - - # Sample List write(imzml_stream, """ """) - - # Software List - OPEN SOURCE write(imzml_stream, """ """) - - # Scan Settings IMAGING write(imzml_stream, """ @@ -506,8 +559,6 @@ function ExportImzml(target_file::String, pixels::Vector{ProcessedPixel}, dims:: """) - - # Instrument Configuration MALDI/TOF write(imzml_stream, """ @@ -525,8 +576,6 @@ function ExportImzml(target_file::String, pixels::Vector{ProcessedPixel}, dims:: """) - - # Data Processing write(imzml_stream, """ @@ -539,18 +588,21 @@ function ExportImzml(target_file::String, pixels::Vector{ProcessedPixel}, dims:: # Run and Spectrum List spectrum_offsets = UInt64[] write(imzml_stream, """ - + """) - # Write each spectrum - for (i, pixel) in enumerate(pixels) - meta = binary_meta[i] - x, y = pixel.coords + # Write each spectrum's metadata + for (i, meta) in enumerate(binary_meta) + x, y = coords[i] spectrum_start = position(imzml_stream) push!(spectrum_offsets, spectrum_start) - write(imzml_stream, """ + # Calculate number of points from byte length + mz_points = meta.mz_length ÷ sizeof(Float64) + int_points = meta.int_length ÷ sizeof(Float32) + + write(imzml_stream, """ @@ -563,14 +615,14 @@ function ExportImzml(target_file::String, pixels::Vector{ProcessedPixel}, dims:: - + - + @@ -603,10 +655,10 @@ function ExportImzml(target_file::String, pixels::Vector{ProcessedPixel}, dims:: return true catch e - @error "Failed to export imzML files" exception=(e, catch_backtrace()) - # Clean up partial files - isfile(ibd_file) && rm(ibd_file, force=true) + @error "Failed to export imzML metadata file" exception=(e, catch_backtrace()) + # Clean up partial .imzML file isfile(target_file) && rm(target_file, force=true) + # Do not delete the .ibd file as it might be useful for debugging return false end end @@ -630,17 +682,15 @@ function ImportMzmlFile(source_file::String, sync_file::String, target_file::Str println("Step 2: Matching acquisition times...") timing_matrix = MatchAcquireTime(sync_file, scans; img_width=img_width, img_height=img_height) - println("Step 3: Converting spectra...") - processed_pixels, (width, height) = ConvertMzmlToImzml(source_file, timing_matrix, scans) + println("Step 3: Converting spectra and writing .ibd file...") + ibd_file = replace(target_file, r"\.imzML$"i => ".ibd") + binary_meta, coords, (width, height) = ConvertMzmlToImzml(source_file, ibd_file, timing_matrix, scans) # Flip image vertically to match R script output - for i in eachindex(processed_pixels) - x, y = processed_pixels[i].coords - processed_pixels[i] = ProcessedPixel((x, height - y + 1), processed_pixels[i].mz, processed_pixels[i].intensity) - end + flipped_coords = [(x, height - y + 1) for (x, y) in coords] - println("Step 4: Exporting to .imzML/.ibd format...") - success = ExportImzml(target_file, processed_pixels, (width, height)) + println("Step 4: Exporting .imzML metadata file...") + success = ExportImzml(target_file, binary_meta, flipped_coords, (width, height)) if success println("Conversion successful: $target_file") diff --git a/src/imzML.jl b/src/imzML.jl index 8b7ec9f..0426109 100644 --- a/src/imzML.jl +++ b/src/imzML.jl @@ -1,4 +1,4 @@ -using Images, Statistics, CairoMakie, DataFrames, Printf, ColorSchemes +using Images, Statistics, CairoMakie, DataFrames, Printf, ColorSchemes, StatsBase # --- Extracted from imzML.jl --- @@ -347,6 +347,38 @@ function load_slices(folder, masses, tolerance) end +""" + get_mz_slice(data::MSIData, mass::Real, tolerance::Real) + +Extracts an image slice for a given m/z value without plotting. +This is a performant function that iterates through spectra once. + +# Returns +- A `Matrix{Float64}` representing the intensity slice. +""" +function get_mz_slice(data::MSIData, mass::Real, tolerance::Real) + width, height = data.image_dims + slice_matrix = zeros(Float64, height, width) # Note: height, width for (y,x) indexing + + _iterate_spectra_fast(data) do spec_idx, mz_array, intensity_array + meta = data.spectra_metadata[spec_idx] + + intensity = find_mass(mz_array, intensity_array, mass, tolerance) + + if intensity > 0.0 + # Populate the matrix using (y, x) indexing + if 1 <= meta.x <= width && 1 <= meta.y <= height + slice_matrix[meta.y, meta.x] = intensity + end + end + end + + # The original function had a NaN replacement, which is good practice to keep. + replace!(slice_matrix, NaN => 0.0) + return slice_matrix +end + + """ plot_slice(msi_data::MSIData, mass::Float64, tolerance::Float64, output_dir::String; stage_name="slice", bins=256) @@ -486,6 +518,30 @@ function set_pixel_depth(img, bounds, depth) return result end +""" + TrIQ(pixMap, depth, prob=0.98) + +Applies TrIQ (Treshold Intensity Quantization) normalization to an image. +This function first computes dynamic intensity range bounds by identifying outliers +based on a cumulative probability, then sets the pixel depth (quantizes intensities) +within these bounds. + +# Arguments +- `pixMap`: The input image matrix (e.g., a slice from `GetMzSliceJl`). +- `depth`: The number of grey levels (bins) to quantize the intensities into. +- `prob`: The target cumulative probability (e.g., 0.98) used to determine outlier thresholds. + +# Returns +- A new image matrix with intensities quantized to the specified depth within the TrIQ bounds. +""" +function TrIQ(pixMap, depth, prob=0.98) + # Compute new dynamic range + bounds = get_outlier_thres(pixMap, prob) + + # Set intensity dynamic range + return set_pixel_depth(pixMap, bounds, depth) +end + """ norm_slices_hist(slices, bins; prob=0.98) @@ -515,11 +571,73 @@ function norm_slices_hist(slices, bins; prob=0.98) return (norm_img=norm_img, bounds=mass_bounds) # bounds is now a vector, one per mass end +""" + quantize_intensity(slice::AbstractMatrix{<:Real}, levels::Integer=256) + +Linearly scales the intensity values in a slice to a specified number of levels. +The output is an array of `UInt8` values. +This is a modernized version of `IntQuantCl`. +""" +function quantize_intensity(slice::AbstractMatrix{<:Real}, levels::Integer=256) + max_val = maximum(slice) + if max_val <= 0 + return zeros(UInt8, size(slice)) + end + + # Scale relative to the absolute maximum value to preserve the zero point. + # The original logic used 'colorLevel' which was the max value (e.g., 255). + scale = (levels - 1) / max_val + + # round is equivalent to floor(x+0.5) for positive numbers. + # clamp is used for robustness against floating point inaccuracies. + image = round.(UInt8, clamp.(slice .* scale, 0, levels - 1)) + + return image +end + function median_filter(img) # 3x3 median filter implementation return mapwindow(median, img, (3, 3)) end + + +""" + downsample_spectrum(mz, intensity, n_points=2000) + +Reduces the number of points in a spectrum for faster plotting, while preserving peaks. +It divides the m/z range into `n_points` bins and keeps only the most intense point from each bin. +""" +function downsample_spectrum(mz::AbstractVector, intensity::AbstractVector, n_points::Integer=2000) + if isempty(mz) || length(mz) <= n_points + return mz, intensity + end + + mz_min, mz_max = extrema(mz) + bin_width = (mz_max - mz_min) / n_points + + # We use a vector of tuples to store the max intensity and its corresponding mz for each bin + # (max_intensity, mz_value) + bins = fill((0.0, 0.0), n_points) + + for i in eachindex(mz) + # Determine the bin for the current point + # Bin indices are 1-based + bin_index = min(n_points, floor(Int, (mz[i] - mz_min) / bin_width) + 1) + + # If the current point's intensity is higher than what's in the bin, replace it + if intensity[i] > bins[bin_index][1] + bins[bin_index] = (intensity[i], mz[i]) + end + end + + # Filter out empty bins and separate the mz and intensity values + final_mz = [b[2] for b in bins if b[1] > 0.0] + final_intensity = [b[1] for b in bins if b[1] > 0.0] + + return final_mz, final_intensity +end + # ============================================================================ # # @@ -678,3 +796,97 @@ function plot_slices(slices, names, masses, output_dir; stage_name, bins=256, dp end return fig end + +""" + save_bitmap(name::String, pixMap::Matrix{UInt8}, colorTable::Vector{UInt32}) + +Saves an 8-bit indexed image as a BMP file. +This is a modernized version of `SaveBitmapCl`. +""" +function save_bitmap(name::String, pixMap::Matrix{UInt8}, colorTable::Vector{UInt32}) + # Get image dimensions + height, width = size(pixMap) + + # Normalize pixel values to stretch contrast, as in the original SaveBitmapCl + minVal, maxVal = extrema(pixMap) + if maxVal > minVal + pixMap = round.(UInt8, 255 * (pixMap .- minVal) ./ (maxVal - minVal)) + end + + # Compute row padding (each row must be a multiple of 4 bytes) + padding = (4 - (width % 4)) % 4 + + # BMP color table must have 256 entries for 8-bit images + fullColorTable = Vector{UInt32}(undef, 256) + if length(colorTable) <= 256 + fullColorTable[1:length(colorTable)] .= colorTable + fullColorTable[length(colorTable)+1:end] .= 0 + else + fullColorTable .= colorTable[1:256] + end + + # Compute file dimensions + offset = 14 + 40 + (256 * 4) # 14(file) + 40(info) + 1024(palette) = 1078 + imgBytes = height * (width + padding) + fileSize = offset + imgBytes + + open(name, "w") do stream + # === File Header (14 bytes) === + write(stream, UInt16(0x4D42)) # "BM" + write(stream, UInt32(fileSize)) + write(stream, UInt16(0)) # Reserved + write(stream, UInt16(0)) # Reserved + write(stream, UInt32(offset)) + + # === Info Header (40 bytes) === + write(stream, UInt32(40)) # Info header size + write(stream, Int32(width)) + write(stream, Int32(height)) # Positive for bottom-up storage + write(stream, UInt16(1)) # Number of color planes + write(stream, UInt16(8)) # Bits per pixel + write(stream, UInt32(0)) # Compression (BI_RGB) + write(stream, UInt32(imgBytes)) + write(stream, Int32(2835)) # Pels per meter X (~72 DPI) + write(stream, Int32(2835)) # Pels per meter Y (~72 DPI) + write(stream, UInt32(256)) # Colors in color table + write(stream, UInt32(0)) # Important colors (0 = all) + + # === Color Table === + write(stream, fullColorTable) + + # === Image Pixels (written bottom-up) === + row_buffer = Vector{UInt8}(undef, width + padding) + row_buffer[width+1:end] .= 0 # pre-fill padding bytes + + for i in height:-1:1 + row_data = @view pixMap[i, :] + row_buffer[1:width] = row_data + write(stream, row_buffer) + end + end +end + +# ******************************************************************** +# Viridis color palette (256 colors) +# ******************************************************************** + +""" + generate_palette(colorscheme, n_colors=256) + +Generates a BMP-compatible UInt32 color palette from a ColorScheme. +""" +function generate_palette(colorscheme, n_colors=256) + palette = Vector{UInt32}(undef, n_colors) + colors = get(colorscheme, range(0, 1, length=n_colors)) + for i in 1:n_colors + c = colors[i] + r = round(UInt8, c.r * 255) + g = round(UInt8, c.g * 255) + b = round(UInt8, c.b * 255) + # BMP color table format is 0x00RRGGBB, written little-endian becomes BB GG RR 00 + palette[i] = (UInt32(r) << 16) | (UInt32(g) << 8) | UInt32(b) + end + return palette +end + +const ViridisPalette = generate_palette(ColorSchemes.viridis) diff --git a/test/readme.md b/test/readme.md index bea2736..f80ee72 100644 --- a/test/readme.md +++ b/test/readme.md @@ -1,21 +1,78 @@ -to run the tests: +# JuliaMSI Test Suite -time julia --project=. test/run_tests.jl +This document provides instructions on how to set up and run the test suite for the `JuliaMSI` package. The tests validate the core functionality of the data processing workflows, including loading, converting, and visualizing mass spectrometry data. -to enable the different test cases, there are boolean variables +## Local Installation and Setup -test1 = true -test2 = true -test3 = true +### 1. Prerequisites +Make sure you have Julia installed (at least version 1.6). If not, we recommend using `juliaup` for installation. +- **juliaup**: [https://github.com/JuliaLang/juliaup](https://github.com/JuliaLang/juliaup) +- **Official binaries**: [https://julialang.org/downloads/](https://julialang.org/downloads/) -set them to true to perform that specific test. +### 2. Download the Project +Download and decompress the repository. You can get it from the official Codeberg page or clone it using Git. +- **Download ZIP**: [https://codeberg.org/LabABI/JuliaMSI/archive/main.zip](https://codeberg.org/LabABI/JuliaMSI/archive/main.zip) +- **Git URL**: `https://codeberg.org/LabABI/JuliaMSI.git` -test 1 consists in mzml validation and spectrum plotting +Ensure you know the directory where the `JuliaMSI` files are located. -test 2 consists in imzml conversion and validation +### 3. Download Test Data +The test script requires example mass spectrometry data files. You can download the required test data from Zenodo: +- **Example Data**: [https://doi.org/10.5281/zenodo.10084132](https://doi.org/10.5281/zenodo.10084132) -test 3 consists in imzml validation and spectrum and slice plotting +Download the data and place it in a known location on your computer. -use the global variables: -TEST_MZML_FILE, SPECTRUM_TO_PLOT, CONVERSION_SOURCE_MZML, CONVERSION_SYNC_FILE, CONVERSION_TARGET_IMZML, TEST_IMZML_FILE, MZ_VALUE_FOR_SLICE, MZ_TOLERANCE, COORDS_TO_PLOT, RESULTS_DIR -to modify input that the test script will recieve. +### 4. Configure the Test Script +Before running the tests, you must edit the `test/run_tests.jl` file to point to the test data you downloaded. + +1. Open `test/run_tests.jl` in a text editor. +2. Locate the `CONFIG: PLEASE FILL IN YOUR FILE PATHS HERE` section. +3. Update the constant variables (e.g., `TEST_MZML_FILE`, `CONVERSION_SOURCE_MZML`, `CONVERSION_SYNC_FILE`) with the absolute paths to the corresponding files on your system. + +## Running the Tests + +1. **Navigate to the Project Directory**: + Open your terminal and use the `cd` command to navigate to the root folder of the `JuliaMSI` project. + ```bash + cd /path/to/your/JuliaMSI + ``` + +2. **Execute the Test Script**: + 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 + ``` + +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`. + +### Enabling and Disabling Test Cases +You can run or skip specific test cases by setting the corresponding boolean variables to `true` or `false`. + +```julia +test1 = true # Runs Test Case 1 +test2 = true # Runs Test Case 2 +test3 = true # Runs Test Case 3 +``` + +- **Test Case 1**: Validates a standard `.mzML` file and generates plots for a single spectrum, the total spectrum, and the average spectrum. +- **Test Case 2**: Tests the conversion of a `.mzML` file (and its corresponding `.txt` sync file) into the `.imzML` format. It then validates the newly created file. +- **Test Case 3**: Validates an existing `.imzML` file. It generates plots for a spectrum at specific coordinates, the total spectrum, the average spectrum, and an ion image slice. + +### Input Configuration Variables +All configuration variables are located in the `CONFIG` section of `test/run_tests.jl`. + +- `TEST_MZML_FILE`: Path to the standard `.mzML` file for Test Case 1. +- `SPECTRUM_TO_PLOT`: The index of the spectrum to plot from the `.mzML` file in Test Case 1. +- `CONVERSION_SOURCE_MZML`: Path to the imaging `.mzML` file to be converted in Test Case 2. +- `CONVERSION_SYNC_FILE`: Path to the `.txt` synchronization file corresponding to `CONVERSION_SOURCE_MZML`. +- `CONVERSION_TARGET_IMZML`: The output path for the `.imzML` file generated in Test Case 2. This is also used as the default input for Test Case 3. +- `TEST_IMZML_FILE`: Path to the `.imzML` file to be tested in Test Case 3. +- `MZ_VALUE_FOR_SLICE`: The m/z value to use for creating an ion image slice in Test Case 3. +- `MZ_TOLERANCE`: The tolerance (+/-) for the `MZ_VALUE_FOR_SLICE` when generating the image. +- `COORDS_TO_PLOT`: An `(X, Y)` tuple specifying the coordinates of the spectrum to plot from the `.imzML` file in Test Case 3. +- `RESULTS_DIR`: The directory where output images will be saved. \ No newline at end of file diff --git a/test/run_tests.jl b/test/run_tests.jl index d96c39f..9b2efb2 100644 --- a/test/run_tests.jl +++ b/test/run_tests.jl @@ -32,7 +32,8 @@ using MSI_src # A regular, non-imaging mzML file. # const TEST_MZML_FILE = "/home/pixel/Documents/Cinvestav_2025/Analisis/mzML/T9_A1.mzML" # const TEST_MZML_FILE = "/home/pixel/Documents/Cinvestav_2025/Analisis/CE4_BF_R1/CE4_BF_R1.mzML" -const TEST_MZML_FILE = "/home/pixel/Documents/Cinvestav_2025/Analisis/CE4_BF_R1/CE4_BF_R1.mzML" +# const TEST_MZML_FILE = "/home/pixel/Documents/Cinvestav_2025/Analisis/CE4_BF_R1/CE4_BF_R1.mzML" +const TEST_MZML_FILE = "/home/pixel/Documents/Cinvestav_2025/Analisis/Imaging_paper_spray/Imaging_paper_spray.mzML" const SPECTRUM_TO_PLOT = 1 # Which spectrum to plot from the file # --- Test Case 2: .mzML + Sync File for Conversion --- @@ -49,15 +50,18 @@ const CONVERSION_TARGET_IMZML = "test/results/converted_mzml.imzML" # --- Test Case 3: Standard .imzML file --- # An existing imzML file (can be the one generated from Case 2). -const TEST_IMZML_FILE = CONVERSION_TARGET_IMZML # The output from case 2 +# const TEST_IMZML_FILE = CONVERSION_TARGET_IMZML # The output from case 2 # const TEST_IMZML_FILE = "/home/pixel/Documents/Cinvestav_2025/Analisis/imzML_AP_SMALDI/HR2MSImouseurinarybladderS096.imzML" # const TEST_IMZML_FILE = "/home/pixel/Documents/Cinvestav_2025/Analisis/Imaging_paper_spray/Imaging_paper_spray.imzML" +const TEST_IMZML_FILE = "/home/pixel/Documents/Cinvestav_2025/Analisis/Imaging prueba Roya 1/royaimg.imzML" # The m/z value to use for creating an image slice. -const MZ_VALUE_FOR_SLICE = 309.06 # BF +# const MZ_VALUE_FOR_SLICE = 309.06 # BF # const MZ_VALUE_FOR_SLICE = 896.0 # HR2MSI # const MZ_VALUE_FOR_SLICE = 76.03 # I PS +const MZ_VALUE_FOR_SLICE = 473 # ROYA # const MZ_TOLERANCE = 0.1 -const MZ_TOLERANCE = 1 +# const MZ_TOLERANCE = 1 +const MZ_TOLERANCE = 0.5 # Coordinates to plot a specific spectrum from imzML const COORDS_TO_PLOT = (50, 50) # Example coordinates (X, Y) @@ -65,8 +69,8 @@ const COORDS_TO_PLOT = (50, 50) # Example coordinates (X, Y) # --- Output Directory --- const RESULTS_DIR = "test/results" -test1 = true -test2 = true +test1 = false +test2 = false test3 = true # =================================================================== @@ -216,8 +220,7 @@ function run_test() end # --- Test Case 3: Process an existing .imzML file --- - println(" -" * "="^20 * " Test Case 3: Processing .imzML " * "="^20) + println("" * "="^20 * " Test Case 3: Processing .imzML " * "="^20) if test3 == true # Run new, stronger validation for imzML validate_msi_data(TEST_IMZML_FILE) @@ -226,6 +229,7 @@ function run_test() if isfile(TEST_IMZML_FILE) # Add spectrum plotting for imzML to match Test Case 1 try + """ # Get the msi data from the imzml println("Plotting a sample spectrum from $TEST_IMZML_FILE...") msi_data = OpenMSIData(TEST_IMZML_FILE) @@ -271,6 +275,8 @@ function run_test() output_path = joinpath(RESULTS_DIR, "test_imzml_average_spectrum.png") save(output_path, fig) println("SUCCESS: Total spectrum plot saved to $output_path") + """ + println("No spectrums tested on this try.") catch e println("ERROR during spectrum plotting in Test Case 3: $e") end @@ -279,7 +285,7 @@ function run_test() try println("Testing plot_slice function on $TEST_IMZML_FILE...") msi_data = OpenMSIData(TEST_IMZML_FILE) - plot_slice(msi_data, MZ_VALUE_FOR_SLICE, MZ_TOLERANCE, RESULTS_DIR, stage_name="test_imzml_single_slice") + @time plot_slice(msi_data, MZ_VALUE_FOR_SLICE, MZ_TOLERANCE, RESULTS_DIR, stage_name="test_imzml_single_slice") # The success message is now inside plot_slice catch e println("ERROR during plot_slice test in Test Case 3: $e")