Implemented imzml and mzml lazy loading, imzml to load spectra matrix, working on mzml to imzml

.gitignore

@@ -4,3 +4,4 @@ public/css/imgOver.png
 !public/css/autogenerated.css
 !public/css/LABI_logo.png
 log/*
+R_original_scripts/

app.jl

@@ -6,7 +6,7 @@ using Libz
 using PlotlyBase
 using CairoMakie
 using Colors
-using julia_mzML_imzML
+# using julia_mzML_imzML
 using Statistics
 using NaturalSort
 using Images
@@ -15,6 +15,7 @@ using NativeFileDialog # Opens the file explorer depending on the OS
 using StipplePlotly
 using Base.Filesystem: mv # To rename files in the system
 include("./julia_imzML_visual.jl")
+include("./src/imzML.jl")
 @genietools
 
 # == Reactive code ==
@@ -151,7 +152,7 @@ include("./julia_imzML_visual.jl")
         ),
         margin=attr(l=0,r=0,t=0,b=0,pad=0)
     )
-    traceImg=PlotlyBase.heatmap(x=[], y=[])
+    traceImg=PlotlyBase.heatmap(x=Vector{Float64}(), y=Vector{Float64}())
     @out plotdataImg=[traceImg]
     @out plotlayoutImg=layoutImg
     # For the image in the comparative view
@@ -180,14 +181,14 @@ include("./julia_imzML_visual.jl")
         margin=attr(l=0,r=0,t=120,b=0,pad=0)
     )
     # Dummy 2D scatter plot
-    traceSpectra=PlotlyBase.scatter(x=[], y=[], mode="lines")
+    traceSpectra=PlotlyBase.stem(x=Vector{Float64}(), y=Vector{Float64}(),marker=attr(size=1, color="blue", opacity=0.1))
     # Create conection to frontend
     @out plotdata=[traceSpectra]
     @out plotlayout=layoutSpectra
     @in xCoord=0
     @in yCoord=0
-    @out xSpectraMz=Float64[]
+    @out xSpectraMz = Vector{Float64}()
-    @out ySpectraMz=Float64[]
+    @out ySpectraMz = Vector{Float64}()
 
     # Interactive plot reactions
     @in data_click=Dict{String,Any}()
@@ -207,7 +208,7 @@ include("./julia_imzML_visual.jl")
         margin=attr(l=0,r=0,t=100,b=0,pad=0)
     )
     # Dummy 2D surface plot
-    traceContour=PlotlyBase.scatter(x=[], y=[], mode="lines")
+    traceContour=PlotlyBase.scatter(x=Vector{Float64}(), y=Vector{Float64}(), mode="lines")
     # Create conection to frontend
     @out plotdataC=[traceContour]
     @out plotlayoutC=layoutContour
@@ -233,7 +234,7 @@ include("./julia_imzML_visual.jl")
     x=1:10
     y=1:10
     z=[sin(i * j / 10) for i in x, j in y]
-    trace3D=PlotlyBase.surface(x=[], y=[], z=[],
+    trace3D=PlotlyBase.surface(x=Vector{Float64}(), y=Vector{Float64}(), z=Matrix{Float64}(undef, 0, 0),
                                 contours_z=attr(
                                     show=true,
                                     usecolormap=true,
@@ -362,7 +363,7 @@ include("./julia_imzML_visual.jl")
         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."
             try
-                spectra=LoadImzml(full_route)
+                spectra=load_imzml(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
@@ -588,7 +589,7 @@ include("./julia_imzML_visual.jl")
             plotdataImg, plotlayoutImg, imgWidth, imgHeight=loadImgPlot(imgInt)
             btnOpticalDisable=false
         else
-            traceImg=PlotlyBase.heatmap(x=[], y=[])
+            traceImg=PlotlyBase.heatmap(x=Vector{Float64}(), y=Vector{Float64}())
             plotdataImg=[traceImg]
             msgimg=""
         end
@@ -615,7 +616,7 @@ include("./julia_imzML_visual.jl")
             plotdataImg, plotlayoutImg, imgWidth, imgHeight=loadImgPlot(imgInt)
             btnOpticalDisable=false
         else
-            traceImg=PlotlyBase.heatmap(x=[], y=[])
+            traceImg=PlotlyBase.heatmap(x=Vector{Float64}(), y=Vector{Float64}())
             plotdataImg=[traceImg]
             msgimg=""
         end
@@ -643,7 +644,7 @@ include("./julia_imzML_visual.jl")
             plotdataImgT, plotlayoutImgT, imgWidth, imgHeight=loadImgPlot(imgIntT)
             btnOpticalDisable=false
         else
-            traceImg=PlotlyBase.heatmap(x=[], y=[])
+            traceImg=PlotlyBase.heatmap(x=Vector{Float64}(), y=Vector{Float64}())
             plotdataImgT=[traceImg]
             msgtriq=""
         end
@@ -670,7 +671,7 @@ include("./julia_imzML_visual.jl")
             plotdataImgT, plotlayoutImgT, imgWidth, imgHeight=loadImgPlot(imgIntT)
             btnOpticalDisable=false
         else
-            traceImg=PlotlyBase.heatmap(x=[], y=[])
+            traceImg=PlotlyBase.heatmap(x=Vector{Float64}(), y=Vector{Float64}())
             plotdataImgT=[traceImg]
             msgtriq=""
         end
@@ -699,7 +700,7 @@ include("./julia_imzML_visual.jl")
             plotdataImgComp, plotlayoutImgComp, _, _=loadImgPlot(imgIntComp)
             btnOpticalDisable=false
         else
-            traceImg=PlotlyBase.heatmap(x=[], y=[])
+            traceImg=PlotlyBase.heatmap(x=Vector{Float64}(), y=Vector{Float64}())
             plotdataImgComp=[traceImg]
             msgimgComp=""
         end
@@ -727,7 +728,7 @@ include("./julia_imzML_visual.jl")
             plotdataImgComp, plotlayoutImgComp, _, _=loadImgPlot(imgIntComp)
             btnOpticalDisable=false
         else
-            traceImg=PlotlyBase.heatmap(x=[], y=[])
+            traceImg=PlotlyBase.heatmap(x=Vector{Float64}(), y=Vector{Float64}())
             plotdataImgComp=[traceImg]
             msgimgComp=""
         end
@@ -755,7 +756,7 @@ include("./julia_imzML_visual.jl")
             plotdataImgTComp, plotlayoutImgTComp, _, _=loadImgPlot(imgIntTComp)
             btnOpticalDisable=false
         else
-            traceImg=PlotlyBase.heatmap(x=[], y=[])
+            traceImg=PlotlyBase.heatmap(x=Vector{Float64}(), y=Vector{Float64}())
             plotdataImgTComp=[traceImg]
             msgtriqComp=""
         end
@@ -783,7 +784,7 @@ include("./julia_imzML_visual.jl")
             plotdataImgTComp, plotlayoutImgTComp, _, _=loadImgPlot(imgIntTComp)
             btnOpticalDisable=false
         else
-            traceImg=PlotlyBase.heatmap(x=[], y=[])
+            traceImg=PlotlyBase.heatmap(x=Vector{Float64}(), y=Vector{Float64}())
             plotdataImgTComp=[traceImg]
             msgtriqComp=""
         end

config/env/global.jl

@@ -1 +1,2 @@
-ENV["GENIE_ENV"] = "prod"
+ENV["GENIE_ENV"] = "dev"
+#ENV["GENIE_ENV"] = "prod"

julia_imzML_visual.jl

@@ -473,7 +473,7 @@ function meanSpectrumPlot(mzmlRoute::String)
         xSpectraMz=spectraMz[1,1]
         ySpectraMz=spectraMz[2,1]
     end
-    trace=PlotlyBase.scatter(x=xSpectraMz, y=ySpectraMz, mode="lines")
+    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
     return plotdata, plotlayout, xSpectraMz, ySpectraMz

ongoing_tasks.txt

@@ -1,27 +0,0 @@
-PENDING
-    Rmsi & julia coherence with image creation and colorbar
-        colorbar revamp CURRENTLY ONGOING
-    Create function that makes the imzML from mzML files ?
-    Possibility to add multiple imzML to process at once ?
-    Add comparison image (rotate, transform, translate, transparency) CURRENTLY ONGOING
-
-DONE
-    Quicker start for julia GUI
-    Better UI
-        refined style of page
-        improved responsiveness
-        UI functionality added
-    Code cleanse
-        better code readability
-    Better blind search for images
-    3d topology plots for images 
-    2d contour plots for images
-    Image flip upside down
-    Measure for the time it takes creating image or plot
-    Rotate vertical images
-    Easier input for files
-    TrIQ default values 
-    Comparative for two views
-    Even faster initial boot and subsectuential boot
-    Multiple spectra plot types
-    Plot creation for spectra Per pixel of image

scripts/build.jl

@@ -0,0 +1,23 @@
+# scripts/build.jl
+using PackageCompiler
+using Pkg
+
+project_dir = dirname(@__DIR__)
+Pkg.activate(project_dir)
+
+# Get all the direct dependencies from Project.toml
+deps = keys(Pkg.project().dependencies)
+sysimage_path = joinpath(project_dir, "build", "JuliaMSI_sysimage.so")
+precompile_script = joinpath(project_dir, "scripts", "precompile.jl")
+
+@info "Starting system image compilation. This may take several minutes..."
+
+create_sysimage(
+    deps;
+    sysimage_path=sysimage_path,
+    precompile_execution_file=precompile_script
+)
+
+@info "Compilation finished!"
+@info "You can now run the fast-starting application using the following command:"
+@info "julia --project=. -J build/JuliaMSI_sysimage.so start_MSI_GUI.jl"

scripts/precompile.jl

@@ -0,0 +1,45 @@
+# scripts/precompile.jl
+using Pkg
+Pkg.activate(dirname(@__DIR__))
+
+using Genie
+using HTTP
+
+@info "Precompiling application..."
+
+# Load the app. Mmap=false is recommended for compilation
+Genie.loadapp(pwd(); Mmap=false)
+
+@info "Starting server for precompilation."
+# Start the server in the background
+server = up(1481, "127.0.0.1", async=true)
+
+# Give the server a moment to start up
+sleep(20)
+
+base_url = "http://127.0.0.1:1481"
+@info "Hitting routes to precompile..."
+
+try
+    # Make a GET request to the home page
+    response = HTTP.get(base_url * "/")
+    @info "GET / -> Status: $(response.status)"
+
+    # === IMPORTANT ===
+    # For best performance, you should add more HTTP requests here
+    # to hit ALL your important routes and API endpoints.
+    # This ensures the code for every page gets compiled.
+    # Example:
+    # HTTP.get(base_url * "/contact")
+    # HTTP.post(base_url * "/submit-data", [], "{\"key\":\"value\"}")
+
+catch e
+    @error "Could not hit routes during precompilation." exception=(e, catch_backtrace())
+
+finally
+    # Stop the server
+    @info "Shutting down server."
+    down()
+end
+
+@info "Precompilation tracing finished."

src/Bitmap.jl

@@ -0,0 +1,399 @@
+# ********************************************************************
+# SaveBitmap
+# ********************************************************************
+"""
+    SaveBitmap( name, pixMap, colorTable )
+
+Save a discretized imzML slice as a bitmap file
+
+# Arguments
+* `      name`: Full path name of target bitmap
+* `    pixMap`: [UInt8] Bidimensional matrix with image info
+* `colorTable`: [UInt32] Vector with RGB colors for each gray level
+
+# Examples
+```julia
+# Saves a bitmap with black and white alternate image squares
+img = zeros(UInt8, 32, 32)
+img[17:32,1:16] .= 1
+img[1:16,17:32] .= 1
+SaveBitmap( "test.bmp", img, UInt32[0, 0x00FFFFFF] )
+```
+
+The following image will be created on your hard disk
+
+![](./Test.bmp)
+"""
+function SaveBitmap( name,
+  pixMap::Array{UInt8,2},
+  colorTable::Array{UInt32,1} )
+
+  # Get image dimensions
+  dim = size( pixMap )
+  if length( dim ) != 2
+    return 0
+  end
+
+  # 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
+
+
+# ********************************************************************
+# FindMass
+#   massVector: mz Vector sorted in ascending order
+#         mass: target mz value
+#    tolerance: bi-axial tolerance to find mass
+# ********************************************************************
+function FindMass( massVector, mass, tolerance )
+
+  index  = Int( 0 )
+  lower  = Int( 1 )
+  higher = Int( length( massVector ) )
+
+  while lower <= higher
+
+    # Compute mid element
+    index = ( lower + higher ) ÷ 2
+
+    # Go to lower portion ?
+    if massVector[ index ] > ( mass + tolerance )
+      higher = index - 1
+      continue
+    end
+
+    # Go to Higher portion?
+    if massVector[ index ] < ( mass - tolerance )
+      lower = index + 1
+      continue
+    end
+
+    # mass was found
+    return index
+  end
+
+  # mass was not found
+  return 0
+end
+
+
+# ********************************************************************
+# GetSlice
+#       imzML: imzML data
+#        mass: mz target value
+#   tolerance: bi-axial tolerance to find mass
+# ********************************************************************
+"""
+    GetSlice( imzML, mass, tolerance )
+
+Extract a mz-image from an imzML image data array loaded with
+ `LoadImzml` function. The resulting image array preserves the same
+ data type of the y-axis vectors in the file.
+
+# Arguments
+* `    imzML`: Image array loaded with `LoadImzml` function
+* `     mass`: mz value for image extraction
+* `tolerance`: Maximum lateral tolerance for a valid mz search
+
+# Examples
+```julia
+# Load DESI MSI Carcinoma  885.55 mz slice
+spectra = LoadImzml( "80TopL, 50TopR, 70BottomL, 60BottomR-centroid.imzML" )
+slice   = GetSlice( spectra, 885.55, 0.005 )
+```
+"""
+function GetSlice( imzML, mass, tolerance )
+
+  # Alloc space for slice
+  width  = maximum( imzML[1,:] )
+  height = maximum( imzML[2,:] )
+  image  = zeros( Float64, width, height )
+
+  for i in 1:size( imzML )[2]
+    index = FindMass( imzML[3,i], mass, tolerance )
+    if index != 0
+     image[ imzML[1,i], imzML[2,i] ] = imzML[4,i][index]
+    end
+  end
+
+  return image
+
+end
+
+
+# ********************************************************************
+# IntQuant, Discretize image amplitude in 0:255 range
+#   slice: Image matriz in measured arbitrary units
+# ********************************************************************
+"""
+    IntQuant( slice )
+
+Zero Memory intensity quantizer. Discretize the image?s continuos
+ intensity range in discrete bins from 0 to 255 gray levels.
+ Returns a UInt8 bidimensional matrix 
+
+# Arguments
+* `slice`: slice returned by the `GetSlice` function
+
+# Examples
+```julia
+# Save DESI MSI Carcinoma 885.55 mz slice as a bitmap file
+spectra = LoadImzml( "80TopL, 50TopR, 70BottomL, 60BottomR-centroid.imzML" ) 
+slice   = GetSlice( spectra, 885.55, 0.005 
+SaveBitmap( "Slice.bmp", IntQuant( slice ), ViridisPalette )
+```
+
+The following image will be created on your hard disk
+
+![](./Slice.bmp)
+"""
+function IntQuant( slice )
+
+  # Compute scale factor for amplitude discretization
+  lower  = minimum( slice )
+  scale  = 255 / 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
+
+
+
+# ********************************************************************
+# TrIQ, Discretize image amplitude in 0:255 range, grouping outliers
+# in the highest bin
+#   slice: Image matrix in measured arbitrary units
+# ********************************************************************
+function GetOutlierThres( slice, prob )
+
+  # Get bin width
+  low = minimum( slice )
+  upp = maximum( slice )
+
+  # Compute histogram's bin count & separation
+  if ( upp - low + 1 ) >= 100
+    bins    = convert( Int32, 100 )
+    step    = ( upp + (upp-low)/(bins-1) - low ) / bins
+  else
+    bins   = convert( Int32, ceil( upp ) - floor( low ) + 1 )
+    step   = 1
+  end
+
+  # Compute histogram
+  histCount = zeros( bins, 1 )
+  for k in slice
+    index = convert( Int64, floor( ( k - low ) / step ) + 1 )
+    histCount[ index ] += 1
+  end
+
+  # Get bin index, that accounts for desired probability
+  prob  = 0.95
+  delta = cumsum( histCount, dims=1 ) / sum( histCount ) .- prob
+  index = findmin( broadcast( abs, delta ) )[2][1]
+
+  # Find max intensity image value
+  key   = low + step * index
+  upp   = -1
+
+  for k in slice[:]
+    if key > k  &&  k > upp
+      upp = k
+    end
+  end
+
+  return [ low, upp ]
+
+end
+
+
+# ********************************************************************
+# Discretize image's intensity using a given intensity range
+#    slice: Image matrix in measured arbitrary units
+#   bounds: max and min intensity discretizing range
+#    depth: number of intensity discrete steps in final image
+# ********************************************************************
+function SetPixelDepth( slice, bounds, depth )
+
+  # Compute intensity bins
+  bins   = depth - 1
+  limits = collect(
+    range( bounds[1], stop = bounds[2], length = depth ) )
+
+  # Reserve memory for output image
+  imgBytes = zeros( UInt8, size( slice ) )
+  nPixels  = length( slice )
+  step     = limits[2] - limits[1]
+
+  # Set intensity depth
+  for k in 1:nPixels
+
+    # Find bin
+    i::Int64 = floor( ( slice[k] - limits[1] ) / step ) + 1
+    if i < bins
+      imgBytes[k] = i + ( slice[k] > limits[i] ) - 1
+    else
+      imgBytes[k] = bins
+    end
+
+  end
+
+  return imgBytes
+end
+
+
+# ********************************************************************
+# Discretize image's intensity removing intensity outliers
+#   slice: Image matrix in measured arbitrary units
+#   depth: Number of intensity discrete steps in final image
+#    prob: Proportion of pixesl to take into account in discretization
+# ********************************************************************
+"""
+    TrIQ( slice, depth, prob = 0.98 )
+
+TrIQ intensity quantizer as described in _DOI 10.7717/peerj-cs.585_.
+ Discretize the image?s continuos intensity range in discrete bins
+ from 0 to 255 gray levels, groupping intensity outiers in the
+ highehst bin. Returns a UInt8 bidimensional matrix 
+
+# Arguments
+* `slice`: slice returned by the `GetSlice` function
+* `depth`: Number ob discrete bins for intensity quantization
+* ` prob`: Cummulative distribute function cutting value
+
+# Examples
+```julia
+# Save DESI MSI Carcinoma 885.55 mz slice as a bitmap file
+spectra = LoadImzml( "80TopL, 50TopR, 70BottomL, 60BottomR-centroid.imzML" ) 
+slice   = GetSlice( spectra, 885.55, 0.005 ) 
+SaveBitmap( "TrIQ.bmp", TrIQ( slice, 256, 0.95 ), ViridisPalette )
+```
+
+The following image will be created on your hard disk
+
+![](./TrIQ.bmp)
+"""
+function TrIQ( slice, depth, prob )
+
+  return  SetPixelDepth(
+    slice,
+    GetOutlierThres( slice, prob ),
+    depth
+  )
+
+end
+
+
+# ********************************************************************
+# Viridis color palette 256 color levels
+# ********************************************************************
+ViridisPalette = UInt32[
+  0x440154, 0x440256, 0x450457, 0x450559,
+  0x46075A, 0x46085C, 0x460A5D, 0x460B5E,
+  0x470D60, 0x470E61, 0x471063, 0x471164,
+  0x471365, 0x481467, 0x481668, 0x481769,
+  0x48186A, 0x481A6C, 0x481B6D, 0x481C6E,
+  0x481D6F, 0x481F70, 0x482071, 0x482173,
+  0x482374, 0x482475, 0x482576, 0x482677,
+  0x482878, 0x482979, 0x472A7A, 0x472C7A,
+  0x472D7B, 0x472E7C, 0x472F7D, 0x46307E,
+  0x46327E, 0x46337F, 0x463480, 0x453581,
+  0x453781, 0x453882, 0x443983, 0x443A83,
+  0x443B84, 0x433D84, 0x433E85, 0x423F85,
+  0x424086, 0x424186, 0x414287, 0x414487,
+  0x404588, 0x404688, 0x3F4788, 0x3F4889,
+  0x3E4989, 0x3E4A89, 0x3E4C8A, 0x3D4D8A,
+  0x3D4E8A, 0x3C4F8A, 0x3C508B, 0x3B518B,
+  0x3B528B, 0x3A538B, 0x3A548C, 0x39558C,
+  0x39568C, 0x38588C, 0x38598C, 0x375A8C,
+  0x375B8D, 0x365C8D, 0x365D8D, 0x355E8D,
+  0x355F8D, 0x34608D, 0x34618D, 0x33628D,
+  0x33638D, 0x32648E, 0x32658E, 0x31668E,
+  0x31678E, 0x31688E, 0x30698E, 0x306A8E,
+  0x2F6B8E, 0x2F6C8E, 0x2E6D8E, 0x2E6E8E,
+  0x2E6F8E, 0x2D708E, 0x2D718E, 0x2C718E,
+  0x2C728E, 0x2C738E, 0x2B748E, 0x2B758E,
+  0x2A768E, 0x2A778E, 0x2A788E, 0x29798E,
+  0x297A8E, 0x297B8E, 0x287C8E, 0x287D8E,
+  0x277E8E, 0x277F8E, 0x27808E, 0x26818E,
+  0x26828E, 0x26828E, 0x25838E, 0x25848E,
+  0x25858E, 0x24868E, 0x24878E, 0x23888E,
+  0x23898E, 0x238A8D, 0x228B8D, 0x228C8D,
+  0x228D8D, 0x218E8D, 0x218F8D, 0x21908D,
+  0x21918C, 0x20928C, 0x20928C, 0x20938C,
+  0x1F948C, 0x1F958B, 0x1F968B, 0x1F978B,
+  0x1F988B, 0x1F998A, 0x1F9A8A, 0x1E9B8A,
+  0x1E9C89, 0x1E9D89, 0x1F9E89, 0x1F9F88,
+  0x1FA088, 0x1FA188, 0x1FA187, 0x1FA287,
+  0x20A386, 0x20A486, 0x21A585, 0x21A685,
+  0x22A785, 0x22A884, 0x23A983, 0x24AA83,
+  0x25AB82, 0x25AC82, 0x26AD81, 0x27AD81,
+  0x28AE80, 0x29AF7F, 0x2AB07F, 0x2CB17E,
+  0x2DB27D, 0x2EB37C, 0x2FB47C, 0x31B57B,
+  0x32B67A, 0x34B679, 0x35B779, 0x37B878,
+  0x38B977, 0x3ABA76, 0x3BBB75, 0x3DBC74,
+  0x3FBC73, 0x40BD72, 0x42BE71, 0x44BF70,
+  0x46C06F, 0x48C16E, 0x4AC16D, 0x4CC26C,
+  0x4EC36B, 0x50C46A, 0x52C569, 0x54C568,
+  0x56C667, 0x58C765, 0x5AC864, 0x5CC863,
+  0x5EC962, 0x60CA60, 0x63CB5F, 0x65CB5E,
+  0x67CC5C, 0x69CD5B, 0x6CCD5A, 0x6ECE58,
+  0x70CF57, 0x73D056, 0x75D054, 0x77D153,
+  0x7AD151, 0x7CD250, 0x7FD34E, 0x81D34D,
+  0x84D44B, 0x86D549, 0x89D548, 0x8BD646,
+  0x8ED645, 0x90D743, 0x93D741, 0x95D840,
+  0x98D83E, 0x9BD93C, 0x9DD93B, 0xA0DA39,
+  0xA2DA37, 0xA5DB36, 0xA8DB34, 0xAADC32,
+  0xADDC30, 0xB0DD2F, 0xB2DD2D, 0xB5DE2B,
+  0xB8DE29, 0xBADE28, 0xBDDF26, 0xC0DF25,
+  0xC2DF23, 0xC5E021, 0xC8E020, 0xCAE11F,
+  0xCDE11D, 0xD0E11C, 0xD2E21B, 0xD5E21A,
+  0xD8E219, 0xDAE319, 0xDDE318, 0xDFE318,
+  0xE2E418, 0xE5E419, 0xE7E419, 0xEAE51A,
+  0xECE51B, 0xEFE51C, 0xF1E51D, 0xF4E61E,
+  0xF6E620, 0xF8E621, 0xFBE723, 0xFDE725 ]
+

src/Common.jl

@@ -0,0 +1,94 @@
+# *******************************************************************
+# VectorConfig, structure with axis decoding instructions
+# *******************************************************************
+mutable struct SpecDim
+  Format::DataType     # Float32, Float64
+  Packed::Bool         # 0: none     1: zlib compression
+  Axis::Int32          # 1: m/z      2: Amplitude
+  Skip::Int64          # Bytes to skip
+end
+
+# *******************************************************************
+# Fill VectorConfig fields from "cvParam" tags
+# *******************************************************************
+function ConfigureSpecDim( stream )
+
+  # Vector configuration
+  axis = SpecDim( Float64, false, 1, 0 )
+
+  # Initial values
+  offset    = position( stream )
+  currLine  = ""
+  matchInfo = RegexMatch
+  
+  while true
+  
+    # Next field
+    currLine  = readline( stream )
+    matchInfo = match( r"^\s*<(cvParam)", currLine )
+    
+    if matchInfo === nothing
+      matchInfo = match( r"^\s*", currLine )
+      axis.Skip = position( stream ) - offset - 
+                  length( currLine ) + length( matchInfo.match )
+      return axis
+    end
+
+    index     = length( matchInfo.captures[1] )
+    matchInfo = GetAttribute( currLine[index:end], "accession" )
+    
+    if matchInfo.captures[1] == "MS:1000515"        # intensity array
+      axis.Axis = 2
+      continue
+      
+    elseif matchInfo.captures[1] == "MS:1000519"    # 32-bit integer
+      axis.Format = Int32
+      continue
+
+    elseif matchInfo.captures[1] == "MS:1000521"    # 32-bit float
+      axis.Format = Float32
+      continue
+      
+    elseif matchInfo.captures[1] == "MS:1000522"    # 64-bit integer
+      axis.Format = Int64
+      continue
+
+    elseif matchInfo.captures[1] == "MS:1000574"		# zlib compresion
+			axis.Packed = true
+			continue
+
+    end
+  end
+end
+
+
+# *******************************************************************
+# Read lines in file until regex match
+#   stream: Source file stream
+#    regex: target tag e.g. "binaryDataArray"
+# *******************************************************************
+function FindTag( stream, regex )
+
+  while true
+    isTag = match( regex, readline( stream ) )
+    if isTag !== nothing
+      return isTag
+    end
+  end  
+end
+
+
+# *******************************************************************
+# Retrieve attribute/value pair
+#   source: String with atr="value" pairs
+#      tag: Optional attribute, e.g. "\\spectrum"
+# *******************************************************************
+function GetAttribute( source, tag = "([^=]+)" )
+
+  # Build regex and retrieve attribute/value pair
+  regStr = Regex( "\\s" * tag * "=\"([^\"]*)\"" )
+  return( match( regStr, source ) )
+  
+end
+
+

src/DataAccess.jl

@@ -0,0 +1,82 @@
+# src/DataAccess.jl
+
+"""
+    load_spectra(fileName::String)
+
+Eagerly loads all spectra from a .mzML or .imzML file into a 2xN matrix.
+This function follows the eager-loading pattern of the original `LoadMzml`
+for backward compatibility and for analyses that require all data in memory.
+
+Each column in the returned matrix represents a single spectrum:
+- Row 1: m/z array
+- Row 2: Intensity array
+
+For .imzML files, the (x,y) spatial coordinates are discarded.
+
+# Arguments
+* `fileName`: Full path to the .mzML or .imzML file.
+
+# Returns
+- An `Array{Any, 2}` containing all spectra.
+"""
+function load_spectra(fileName::String)
+    if endswith(lowercase(fileName), ".mzml")
+        # Detected .mzML file. Using the existing eager loader.
+        return LoadMzml(fileName)
+
+    elseif endswith(lowercase(fileName), ".imzml")
+        # Detected .imzML file. Eagerly loading all spectra from .ibd.
+        
+        # Lazily load metadata to get offsets and counts
+        imzml_data = load_imzml(fileName)
+        
+        try
+            num_spectra = length(imzml_data.spectra_metadata)
+            if num_spectra == 0
+                return Array{Any}(undef, (2, 0))
+            end
+
+            spectra_matrix = Array{Any}(undef, (2, num_spectra))
+            
+            hIbd = imzml_data.ibd_handle
+            mz_format = imzml_data.mz_format
+            int_format = imzml_data.intensity_format
+
+            # For performance, pre-allocate one large buffer for reading.
+            # This avoids re-allocating memory for each spectrum inside the loop.
+            max_points = maximum(meta.mz_count for meta in imzml_data.spectra_metadata)
+            mz_array_buffer = Array{mz_format}(undef, max_points)
+            intensity_array_buffer = Array{int_format}(undef, max_points)
+
+            for i in 1:num_spectra
+                meta = imzml_data.spectra_metadata[i]
+                
+                # Create a view into the buffer with the correct size for this spectrum
+                current_mz_array = view(mz_array_buffer, 1:meta.mz_count)
+                current_int_array = view(intensity_array_buffer, 1:meta.intensity_count)
+
+                # Read binary data directly into the sized views
+                seek(hIbd, meta.mz_offset)
+                read!(hIbd, current_mz_array)
+                
+                seek(hIbd, meta.intensity_offset)
+                read!(hIbd, current_int_array)
+
+                # Store a copy in the final matrix. A copy is necessary because
+                # the buffer is overwritten in the next iteration.
+                spectra_matrix[1, i] = copy(current_mz_array)
+                spectra_matrix[2, i] = copy(current_int_array)
+            end
+            
+            return spectra_matrix
+        finally
+            # Ensure the .ibd file handle is closed, as this is an eager load.
+            if isopen(imzml_data.ibd_handle)
+                close(imzml_data.ibd_handle)
+            end
+        end
+
+    else
+        error("Unsupported file type. Please provide a .mzML or .imzML file.")
+    end
+end

src/MzmlConverter.jl

@@ -0,0 +1,270 @@
+# 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.
+It replicates the functionality of the project's original R scripts that use MALDIquant.
+"""
+
+using CSV, DataFrames
+
+# Note: This file assumes that `LoadSpectra` from DataAccess.jl and the necessary
+# structs/functions from ParserHelpers.jl, mzML.jl, and imzML.jl are available
+# in the execution context (e.g., included in the main app.jl).
+
+# A struct to hold the processed pixel data before exporting
+struct ProcessedPixel
+    coords::Tuple{Int, Int}
+    mz::Vector{Float64} # Assuming m/z is consistent, can be optimized later
+    intensity::Vector{Float32}
+end
+
+"""
+    GetMzmlScanTime(fileName::String)
+
+Parses a .mzML file to extract the scan start time for each spectrum.
+This is a Julia implementation of the `GetMzmlScanTime` function from the R scripts.
+
+# Arguments
+* `fileName`: Path to the .mzML file.
+
+# Returns
+- A `Matrix{Int64}` where each row is `[spectrum_index, time_in_milliseconds]`.
+"""
+function GetMzmlScanTime(fileName::String)
+    spec_count = 0
+    times = Tuple{Int64, Int64}[]
+
+    open(fileName, "r") do stream
+        file_content = read(stream, String)
+        
+        # First, find the total number of spectra to pre-allocate
+        spectrum_list_match = match(r"<spectrumList count=\"(\d+)\"", file_content)
+        if spectrum_list_match === nothing
+            error("Could not find <spectrumList> tag or count attribute.")
+        end
+        total_spectra = parse(Int, spectrum_list_match.captures[1])
+        sizehint!(times, total_spectra)
+
+        # Iterate over each spectrum block
+        for spectrum_match in eachmatch(r"<spectrum index=\"(\d+)\"[^>]*>.*?<\/spectrum>", file_content, overlay=true)
+            spec_block = spectrum_match.match
+            # The R script appears to treat the 0-based mzML index as 1-based, so we add 1.
+            spec_index = parse(Int, spectrum_match.captures[1]) + 1
+
+            # Find the scan start time within the spectrum block
+            time_match = match(r"<cvParam accession=\"MS:1000016\".*?value=\"([\d\.]+)\"", spec_block)
+            
+            if time_match !== nothing
+                # The value is in minutes. Convert to milliseconds and round.
+                time_in_minutes = parse(Float64, time_match.captures[1])
+                time_in_ms = round(Int64, time_in_minutes * 60000)
+                push!(times, (spec_index, time_in_ms))
+            end
+        end
+        
+        # Normalize times so the first scan is at t=0
+        if !isempty(times)
+            sort!(times, by = x -> x[1]) # Ensure scans are sorted by index
+            first_time = times[1][2]
+            for i in eachindex(times)
+                times[i] = (times[i][1], times[i][2] - first_time)
+            end
+        end
+    end # open
+
+    # Convert vector of tuples to a matrix
+    return permutedims(hcat(collect.(times)...))
+end
+
+"""
+    MatchAcquireTime(sync_file_path::String, scans::Matrix{Int64})
+
+Correlates pixel acquisition times from a synchronization file with scan
+acquisition times from an mzML file.
+
+This is a Julia implementation of the `MatchAcquireTime` function from the R scripts.
+
+# Arguments
+* `sync_file_path`: Path to the synchronization file (.txt).
+* `scans`: A matrix of scan times, as returned by `GetMzmlScanTime`.
+
+# Returns
+- A `Matrix{Int64}` where each row is `[pixel_index, pixel_time_ms, first_scan_index, last_scan_index]`.
+"""
+function MatchAcquireTime(sync_file_path::String, scans::Matrix{Int64})
+    if !isfile(sync_file_path)
+        error("Synchronization file not found: $sync_file_path")
+    end
+
+    # Load xyz stage pixel delays from the sync file, skipping the first 2 header lines.
+    pixel_df = CSV.read(sync_file_path, DataFrame, header=false, skipto=3)
+    pixel = Matrix(pixel_df)
+    
+    # Normalize pixel times so the first is at t=0
+    pixel[:, 2] .-= pixel[1, 2]
+
+    num_pixels = size(pixel, 1)
+    num_scans = size(scans, 1)
+    
+    index_matrix = zeros(Int64, num_pixels, 2)
+    
+    scan_row = 1
+    for pixel_row in 1:(num_pixels - 1)
+        start_scan = scan_row
+        
+        # Find the scan that ends after the *next* pixel starts
+        while scan_row <= num_scans && scans[scan_row, 2] < pixel[pixel_row + 1, 2]
+            scan_row += 1
+        end
+        
+        index_matrix[pixel_row, 1] = start_scan
+        index_matrix[pixel_row, 2] = scan_row - 1
+        
+        # The next pixel's scans start from the one that crossed the boundary
+        scan_row = max(1, scan_row - 1)
+    end
+    
+    # Assign scans for the last pixel
+    index_matrix[num_pixels, 1] = scan_row
+    index_matrix[num_pixels, 2] = num_scans
+    
+    # Handle cases where some pixels have no scans by creating an empty range
+    for i in 1:num_pixels
+        if index_matrix[i, 2] < index_matrix[i, 1]
+            index_matrix[i, 2] = index_matrix[i, 1] - 1
+        end
+    end
+
+    return hcat(pixel, index_matrix)
+end
+
+
+function RenderPixel(pixel_info, scans, spectra, scan_time_deltas, pixel_time_deltas)
+    first_scan = pixel_info[3]
+    last_scan = pixel_info[4]
+    num_actions = last_scan - first_scan
+
+    # Use the m/z array of the first scan as the reference
+    mz_array = spectra[1, first_scan]
+    new_intensity = zeros(Float32, length(mz_array))
+
+    if num_actions == 0 # Single scan contributes to the pixel
+        scale = pixel_time_deltas[pixel_info[1]] / scan_time_deltas[first_scan]
+        new_intensity .+= spectra[2, first_scan] .* scale
+
+    elseif num_actions == 1 # Two partial scans contribute
+        # First partial scan
+        scale1 = (scans[last_scan, 2] - pixel_info[2]) / scan_time_deltas[first_scan]
+        new_intensity .+= spectra[2, first_scan] .* scale1
+
+        # Second partial scan
+        next_pixel_time = pixel_info[2] + pixel_time_deltas[pixel_info[1]]
+        scale2 = (next_pixel_time - scans[last_scan, 2]) / scan_time_deltas[last_scan]
+        new_intensity .+= spectra[2, last_scan] .* scale2
+
+    elseif num_actions > 1 # Multiple scans contribute
+        # First partial scan
+        scale1 = (scans[first_scan + 1, 2] - pixel_info[2]) / scan_time_deltas[first_scan]
+        new_intensity .+= spectra[2, first_scan] .* scale1
+
+        # Full scans in the middle
+        for i in (first_scan + 1):(last_scan - 1)
+            new_intensity .+= spectra[2, i]
+        end
+
+        # Last partial scan
+        next_pixel_time = pixel_info[2] + pixel_time_deltas[pixel_info[1]]
+        scale2 = (next_pixel_time - scans[last_scan, 2]) / scan_time_deltas[last_scan]
+        new_intensity .+= spectra[2, last_scan] .* scale2
+    end
+
+    return (mz_array, new_intensity)
+end
+
+function ConvertMzmlToImzml(source_file::String, timing_matrix::Matrix{Int64}, scans::Matrix{Int64})
+    # Get image width from timing info
+    width = findfirst(i -> timing_matrix[i, 1] - timing_matrix[i-1, 1] != 1, 2:size(timing_matrix, 1))
+    height = size(timing_matrix, 1) ÷ width
+
+    # Load the full mzML data
+    spectra = LoadSpectra(source_file)
+
+    # Pre-calculate time deltas
+    scan_time_deltas = diff(scans[:, 2])
+    pixel_time_deltas = diff(timing_matrix[:, 2])
+    # Append a final delta for the last element
+    push!(scan_time_deltas, scan_time_deltas[end])
+    push!(pixel_time_deltas, pixel_time_deltas[end])
+
+    processed_pixels = ProcessedPixel[]
+    sizehint!(processed_pixels, size(timing_matrix, 1))
+
+    for i in 1:size(timing_matrix, 1)
+        pixel_info = timing_matrix[i, :]
+        first_scan = pixel_info[3]
+        last_scan = pixel_info[4]
+
+        # Skip if there are no scans for this pixel
+        (first_scan > last_scan) && continue
+
+        # Calculate coordinates
+        x = ((pixel_info[1] - 1) % width) + 1
+        y = fld(pixel_info[1] - 1, width) + 1
+
+        mz, intensity = RenderPixel(pixel_info, scans, spectra, scan_time_deltas, pixel_time_deltas)
+        push!(processed_pixels, ProcessedPixel((x, y), mz, intensity))
+    end
+
+    return processed_pixels, (width, height)
+end
+
+function ExportImzml(target_file::String, pixels::Vector{ProcessedPixel}, dims::Tuple{Int, Int})
+    @warn "ExportImzml is not yet implemented. The processed pixel data has been generated but not saved to .imzML/.ibd files."
+    # TODO: Implement the logic to write the .imzML (XML) and .ibd (binary) files.
+    # 1. Open target_file.imzML and target_file.ibd for writing.
+    # 2. Write all m/z and intensity arrays sequentially to the .ibd file, tracking offsets and lengths.
+    # 3. Write the .imzML XML structure, including:
+    #    - Boilerplate headers
+    #    - <fileDescription>, <referenceableParamGroupList>, <scanSettingsList>
+    #    - A <spectrumList> with a <spectrum> for each pixel.
+    #    - Each <spectrum> must contain cvParams for x/y coords and external data pointers to the .ibd file.
+    return false
+end
+
+"""
+    ImportMzmlFile(source_file::String, sync_file::String, target_file::String)
+
+Main workflow function to convert a .mzML file to an .imzML file.
+
+# Arguments
+* `source_file`: Path to the input .mzML file.
+* `sync_file`: Path to the synchronization text file.
+* `target_file`: Path for the output .imzML file (the .ibd will be named accordingly).
+"""
+function ImportMzmlFile(source_file::String, sync_file::String, target_file::String)
+    println("Step 1: Getting scan times from .mzML file...")
+    scans = GetMzmlScanTime(source_file)
+
+    println("Step 2: Matching acquisition times...")
+    timing_matrix = MatchAcquireTime(sync_file, scans)
+
+    println("Step 3: Converting spectra...")
+    processed_pixels, (width, height) = ConvertMzmlToImzml(source_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
+
+    println("Step 4: Exporting to .imzML/.ibd format...")
+    success = ExportImzml(target_file, processed_pixels, (width, height))
+
+    if success
+        println("Conversion successful: $target_file")
+    else
+        println("Conversion failed. Exporting is not yet implemented.")
+    end
+    return success
+end

src/ParserHelpers.jl

@@ -0,0 +1,122 @@
+# src/ParserHelpers.jl
+
+"""
+This file provides common helper functions for parsing mzML and imzML files.
+"""
+
+# ============================================================================
+#
+# Data Structures and Low-Level Parser Helpers
+#
+# ============================================================================
+
+"""
+    SpecDim
+
+A struct to hold configuration for a spectral data axis (e.g., m/z or intensity).
+"""
+mutable struct SpecDim
+    Format::Type
+    Packed::Bool
+    Axis::Int
+    Skip::Int
+end
+
+"""
+    find_tag(stream, regex::Regex)
+
+Reads a stream line-by-line until a line matches the provided regex.
+
+# Returns
+- A `RegexMatch` object if a match is found, otherwise throws an error.
+"""
+function find_tag(stream, regex::Regex)
+    while !eof(stream)
+        line = readline(stream)
+        isTag = match(regex, line)
+        if isTag !== nothing
+            return isTag
+        end
+    end
+    error("Tag not found for regex: $regex")
+end
+
+"""
+    get_attribute(source::String, tag::String = "([^=]+)")
+
+Retrieves an attribute's value from an XML tag string.
+
+# Returns
+- A `RegexMatch` object containing the attribute and its value.
+"""
+function get_attribute(source::AbstractString, tag::String = "([^=]+)")
+    # Construct the regex pattern string
+    pattern_str = "\s" * tag * "=\"([^"]*)\""
+    regStr = Regex(pattern_str)
+    return match(regStr, source)
+end
+
+"""
+    configure_spec_dim(stream)
+
+Fills a `SpecDim` struct by parsing `cvParam` tags from the stream.
+"""
+function configure_spec_dim(stream)
+    axis = SpecDim(Float64, false, 1, 0)
+    offset = position(stream)
+
+    while !eof(stream)
+        currLine = readline(stream)
+        matchInfo = match(r"^\s*<(cvParam)", currLine)
+
+        if matchInfo === nothing
+            matchInfo = match(r"^\s*", currLine)
+            axis.Skip = position(stream) - offset - length(currLine) + length(matchInfo.match)
+            return axis
+        end
+
+        index = length(matchInfo.captures[1])
+        attr_match = get_attribute(currLine[index:end], "accession")
+
+        if attr_match !== nothing
+            accession = attr_match.captures[1]
+            if accession == "MS:1000515"       # intensity array
+                axis.Axis = 2
+            elseif accession == "MS:1000519"   # 32-bit integer
+                axis.Format = Int32
+            elseif accession == "MS:1000521"   # 32-bit float
+                axis.Format = Float32
+            elseif accession == "MS:1000522"   # 64-bit integer
+                axis.Format = Int64
+            elseif accession == "MS:1000574"   # zlib compression
+                axis.Packed = true
+            end
+        end
+    end
+    return axis # Should be unreachable if file is well-formed
+end
+
+# ============================================================================
+#
+# Data Structures and Helpers for mzML lazy parsing
+#
+# ============================================================================
+
+mutable struct CVParams
+    format::Type
+    is_compressed::Bool
+    axis_type::Symbol
+end
+
+function update_cv_params!(params::CVParams, acc::String)
+    if acc == "MS:1000514"; params.axis_type = :mz;
+    elseif acc == "MS:1000515"; params.axis_type = :intensity;
+    elseif acc == "MS:1000519"; params.format = Int32;
+    elseif acc == "MS:1000521"; params.format = Float32;
+    elseif acc == "MS:1000522"; params.format = Int64;
+    elseif acc == "MS:1000523"; params.format = Float64;
+    elseif acc == "MS:1000574"; params.is_compressed = true;
+    end
+end
+
+

src/imzML.jl

@@ -0,0 +1,699 @@
+
+using Images, ImageFiltering, StatsBase, Statistics, CairoMakie, ColorSchemes, DataFrames, CSV, Printf, Dates
+include("ParserHelpers.jl")
+
+# --- Lazy Loading Data Structures ---
+
+# Struct to hold metadata for a single spectrum, for lazy loading.
+struct SpectrumMetadata
+    x::Int32
+    y::Int32
+    mz_offset::Int64
+    intensity_offset::Int64
+    mz_count::Int32
+    intensity_count::Int32
+end
+
+# Struct to hold the parsed imzML metadata and the handle to the .ibd file.
+mutable struct ImzMLData
+    ibd_handle::IO
+    mz_format::Type
+    intensity_format::Type
+    spectra_metadata::Vector{SpectrumMetadata}
+    width::Int
+    height::Int
+    mz_is_first::Bool # To know the order in the ibd file
+
+    function ImzMLData(ibd_handle::IO, mz_format::Type, intensity_format::Type, spectra_metadata::Vector{SpectrumMetadata}, width::Int, height::Int, mz_is_first::Bool)
+        obj = new(ibd_handle, mz_format, intensity_format, spectra_metadata, width, height, mz_is_first)
+        finalizer(obj) do o
+            if isopen(o.ibd_handle)
+                close(o.ibd_handle)
+            end
+        end
+        return obj
+    end
+end
+
+
+# --- Extracted from imzML.jl ---
+
+"""
+This file provides a library for parsing `.imzML` and `.ibd` files in pure Julia.
+It is intended to be included by a parent script.
+
+Core Functions:
+- `load_imzml`: The main function that orchestrates the parsing.
+- Helper functions for reading XML metadata and binary spectral data.
+"""
+
+
+# ============================================================================
+#
+#
+# imzML Parser Implementation
+#
+# ============================================================================
+
+"""
+    axes_config_img(stream)
+
+Determines the storage order of the m/z and intensity arrays.
+"""
+function axes_config_img(stream)
+    tag = find_tag(stream, r"^\s*<(referenceableParamGroup )")
+    value = get_attribute(tag.captures[1], "intensityArray")
+    order = 1 + (value !== nothing)
+
+    axis = Array{SpecDim,1}(undef, 2)
+    axis[order] = configure_spec_dim(stream)
+
+    find_tag(stream, r"^\s*<(referenceableParamGroup )")
+    axis[xor(order, 3)] = configure_spec_dim(stream)
+    return axis
+end
+
+"""
+    get_img_dimensions(stream)
+
+Reads the maximum X and Y dimensions and total spectrum count from the metadata.
+"""
+function get_img_dimensions(stream)
+    find_tag(stream, r"^\s*<(scanSettings )")
+    n = 2
+    dim = [0, 0, 0]
+
+    while n > 0 && !eof(stream)
+        currLine = readline(stream)
+        if occursin("<cvParam", currLine)
+            accession_match = get_attribute(currLine, "accession")
+            if accession_match !== nothing
+                accession = accession_match.captures[1]
+                if accession == "IMS:1000042" || accession == "IMS:1000043"
+                    value_match = get_attribute(currLine, "value")
+                    if value_match !== nothing
+                        axis_idx = (accession == "IMS:1000042") ? 1 : 2
+                        dim[axis_idx] = parse(Int32, value_match.captures[1])
+                        n -= 1
+                    end
+                end
+            end
+        end
+    end
+
+    count_tag = find_tag(stream, r"^\s*<spectrumList(.+)")
+    count_match = get_attribute(count_tag.captures[1], "count")
+    dim[3] = parse(Int32, count_match.captures[1])
+    return dim
+end
+
+"""
+    get_spectrum_tag_offset(stream)
+
+Calculates the character offset within a `<spectrum>` tag, ignoring attribute values.
+"""
+function get_spectrum_tag_offset(stream)
+    offset = position(stream)
+    tag = find_tag(stream, r"^\s*<spectrum (.+)")
+    first = 1
+
+    while true
+        value = match(r"[^=]+\"([^\"]+)\"", tag.captures[1][first:end])
+        if value === nothing
+            break
+        end
+        first += value.offsets[1] + length(value.captures[1])
+        offset += length(value.captures[1])
+    end
+    return offset
+end
+
+"""
+    get_spectrum_attributes(stream, hIbd)
+
+Reads metadata to determine the byte offsets and data types for reading spectra.
+"""
+function get_spectrum_attributes(stream, hIbd)
+    skip = Vector{UInt32}(undef, 8)
+    offset = get_spectrum_tag_offset(stream)
+
+    tag = find_tag(stream, r" accession=\"IMS:100005(\d)\"(.+)")
+    skip[1] = tag.captures[1][1] - '0' + 1
+    skip[2] = xor(skip[1], 3)
+
+    value = match(r" value=\"(\d+)\".+", tag.captures[2])
+    skip[5] = position(stream) - offset - length(value.match) - 2
+    value = match(r" value=\"\d+\".+", readline(stream))
+    skip[6] = value.offset
+
+    offset = position(stream)
+    tag = find_tag(stream, r"^\s*<referenceableParamGroupRef(.+)")
+    value = get_attribute(tag.captures[1], "ref")
+    skip[3] = (value.captures[1] == "intensityArray") + 3
+    skip[4] = xor(skip[3], 7)
+
+    k = 2
+    while k != 0
+        tag = find_tag(stream, r" accession=\"IMS:100010(\d)\"(.+)")
+        accession_val = tag.captures[1][1]
+        if accession_val == '2'
+            value = get_attribute(tag.match, "value")
+            seek(hIbd, parse(Int64, value.captures[1]))
+            k -= 1
+        elseif accession_val == '3'
+            skip[7] = position(stream) - offset - length(tag.match)
+            offset = position(stream)
+            k -= 1
+        end
+    end
+
+    find_tag(stream, r"^\s*</spectrum>")
+    skip[8] = position(stream) - offset
+    return skip
+end
+
+
+"""
+    load_imzml(file_path::String)
+Main function to parse an `.imzML`/.ibd file pair and prepare for lazy loading.
+It parses the metadata from the `.imzML` file, including spectrum coordinates and
+byte offsets for data in the `.ibd` file, but does not load the spectral data itself.
+
+# Returns
+- An `ImzMLData` struct containing the open `.ibd` file handle and all necessary metadata
+  for on-demand data loading.
+"""
+function load_imzml(file_path::String)
+    if !isfile(file_path)
+        error("Provided path is not a file: $(file_path)")
+    end
+
+    ibd_path = replace(file_path, ".imzML" => ".ibd")
+    if !isfile(ibd_path)
+        error("Corresponding .ibd file not found for: $(file_path)")
+    end
+
+    stream = open(file_path)
+    hIbd = open(ibd_path)
+
+    # We don't use a try/finally block because we need to return the open hIbd handle.
+    # The finalizer on the ImzMLData struct will be responsible for closing it.
+
+    axis = axes_config_img(stream)
+    imgDim = get_img_dimensions(stream)
+    width, height, num_spectra = imgDim
+
+    # Determine which format is for m/z and which for intensity from the referenceableParamGroups
+    mz_config_idx = findfirst(a -> a.Axis == 1, axis)
+    int_config_idx = findfirst(a -> a.Axis == 2, axis)
+    mz_format = axis[mz_config_idx].Format
+    intensity_format = axis[int_config_idx].Format
+
+    # The original parsing logic is highly optimized and relies on fixed offsets within the XML structure.
+    # We reuse it to get the initial .ibd file offset and the structural attributes of the <spectrum> tags.
+    start_of_spectra_xml = position(stream)
+    attr = get_spectrum_attributes(stream, hIbd)
+    
+    # get_spectrum_attributes moves the hIbd pointer to the start of the first spectrum's binary data.
+    current_ibd_offset = position(hIbd)
+
+    # We need to re-scan the spectra list from the beginning to get metadata for each spectrum.
+    seek(stream, start_of_spectra_xml)
+
+    spectra_metadata = Vector{SpectrumMetadata}(undef, num_spectra)
+    
+    # attr[3] is 4 if the first binary data array is intensity, 3 if it's m/z.
+    mz_is_first = attr[3] == 3
+
+    for k in 1:num_spectra
+        # These skips are based on the structure of the first spectrum tag, assuming all are identical.
+        skip(stream, attr[5]) # Skip to X coordinate value
+        val_tag = find_tag(stream, r"value=\"(\d+)\"")
+        x = parse(Int32, val_tag.captures[1])
+
+        skip(stream, attr[6]) # Skip to Y coordinate value
+        val_tag = find_tag(stream, r"value=\"(\d+)\"")
+        y = parse(Int32, val_tag.captures[1])
+
+        skip(stream, attr[7]) # Skip to array length value
+        val_tag = find_tag(stream, r"value=\"(\d+)\"")
+        nPoints = parse(Int32, val_tag.captures[1])
+
+        mz_len_bytes = nPoints * sizeof(mz_format)
+        int_len_bytes = nPoints * sizeof(intensity_format)
+
+        local mz_offset, int_offset
+        if mz_is_first
+            mz_offset = current_ibd_offset
+            int_offset = mz_offset + mz_len_bytes
+        else
+            int_offset = current_ibd_offset
+            mz_offset = int_offset + int_len_bytes
+        end
+
+        spectra_metadata[k] = SpectrumMetadata(x, y, mz_offset, int_offset, nPoints, nPoints)
+        
+        # Advance the offset for the next spectrum's data.
+        current_ibd_offset += mz_len_bytes + int_len_bytes
+        
+        skip(stream, attr[8]) # Skip to the end of the spectrum tag
+    end
+    
+    # Close the .imzML file stream, but leave the .ibd stream open for lazy reading.
+    close(stream)
+
+    return ImzMLData(hIbd, mz_format, intensity_format, spectra_metadata, width, height, mz_is_first)
+end
+
+# --- End of content from imzML.jl ---
+
+
+# --- Start of content from Imaging_Normalization.jl ---
+
+# =============================================================================
+#
+# Image Slice Extraction
+#
+# =============================================================================
+
+"""
+    find_mass(mz_array, intensity_array, target_mass, tolerance)
+
+Finds the intensity of the most intense peak within a mass tolerance window.
+This modernized version is more robust than a simple binary search as it
+correctly handles multiple peaks within the tolerance window.
+
+# Returns
+- The intensity (`Float64`) of the peak if found, otherwise `0.0`.
+"""
+function find_mass(mz_array, intensity_array, target_mass, tolerance)
+    lower_bound = target_mass - tolerance
+    upper_bound = target_mass + tolerance
+
+    max_intensity = 0.0
+    found = false
+
+    # Iterate through the spectrum to find the highest intensity peak in the window
+    for i in eachindex(mz_array)
+        if lower_bound <= mz_array[i] <= upper_bound
+            if intensity_array[i] > max_intensity
+                max_intensity = intensity_array[i]
+                found = true
+            end
+        end
+    end
+
+    return found ? max_intensity : 0.0
+end
+
+"""
+    load_slices(folder, masses, tolerance)
+
+Loads image slices for multiple masses from all `.imzML` files in a directory.
+This function is optimized to read the binary data file (`.ibd`) only once per file.
+It iterates through each spectrum, reads its data, and finds all target masses
+before proceeding to the next spectrum, minimizing disk I/O.
+"""
+function load_slices(folder, masses, tolerance)
+    files = filter(f -> endswith(f, ".imzML"), readdir(folder, join=true))
+    if isempty(files)
+        @warn "No .imzML files found in the specified directory: $folder"
+        return (Array{Any}(undef, 0, 0), String[])
+    end
+    n_files = length(files)
+    n_slices = length(masses)
+
+    img_list = Array{Any}(undef, n_files, n_slices)
+    names = String[]
+
+    for (i, file) in enumerate(files)
+        name = basename(file)
+        push!(names, name)
+        @info "Processing $(i)/$(n_files): $(name)"
+        
+        imzML_data = @time load_imzml(file) # This is fast, it's fine.
+
+        # Create empty images for all slices for the current file
+        current_file_slices = [zeros(Float64, imzML_data.width, imzML_data.height) for _ in 1:n_slices]
+
+        hIbd = imzML_data.ibd_handle
+        mz_format = imzML_data.mz_format
+        int_format = imzML_data.intensity_format
+
+        # Pre-allocate arrays to be reused for each spectrum
+        mz_array = Array{mz_format}(undef, 0)
+        intensity_array = Array{int_format}(undef, 0)
+
+        # Iterate through each spectrum ONCE
+        for meta in imzML_data.spectra_metadata
+            # Resize arrays if the number of points in the spectrum has changed
+            if length(mz_array) != meta.mz_count
+                resize!(mz_array, meta.mz_count)
+                resize!(intensity_array, meta.intensity_count)
+            end
+
+            # Read spectrum data ONCE from the .ibd file
+            seek(hIbd, meta.mz_offset)
+            read!(hIbd, mz_array)
+            seek(hIbd, meta.intensity_offset)
+            read!(hIbd, intensity_array)
+
+            # Now, check for all masses of interest in this single spectrum
+            for (j, mass) in enumerate(masses)
+                intensity = find_mass(mz_array, intensity_array, mass, tolerance)
+                if intensity > 0.0
+                    if 1 <= meta.x <= imzML_data.width && 1 <= meta.y <= imzML_data.height
+                        current_file_slices[j][meta.x, meta.y] = intensity
+                    end
+                end
+            end
+        end # end of spectra loop
+
+        # Assign the generated images to the main list
+        for j in 1:n_slices
+            img_list[i, j] = current_file_slices[j]
+        end
+
+    end # end of files loop
+    
+    return (img_list, names)
+end
+
+# ============================================================================
+#
+#
+# Image Processing and Normalization
+#
+# ============================================================================
+
+function get_outlier_thres(img, prob=0.98)
+    # DO NOT filter zeros. Use all pixel values like R does.
+    int_values = vec(img)
+    low = minimum(int_values)
+    upp = maximum(int_values)
+
+    # Create histogram bins
+    if upp - low + 1 >= 100
+        bins = 100
+        brk = range(low, stop=upp + (upp - low)/(bins - 1), length=bins + 1)
+    else
+        brk = collect(floor(low):(ceil(upp) + 1))
+    end
+
+    # Compute histogram with RIGHT-CLOSED = FALSE to mimic R's right=FALSE
+    # This makes intervals [a, b) instead of the default [a, b]
+    h = fit(Histogram, int_values, brk, closed=:left) # Key change: closed=:left
+
+    # Check if histogram is valid
+    if isempty(h.weights) || sum(h.weights) == 0
+        return (low, upp)
+    end
+
+    # Replicate R's algorithm exactly
+    cum_counts = cumsum(h.weights) / sum(h.weights)
+    # Calculate the difference from the target probability
+    top = prob .- cum_counts
+    delta = abs.(top)
+    # Find the index where the difference is minimized
+    min_delta_index = findfirst(x -> x == minimum(delta), delta)
+    
+    # R adds 1 to this index: index <- 1 + which(...)[1]
+    target_bin_index = min_delta_index + 1
+    # Get the upper edge of the target bin
+    target_bin_upper_edge = h.edges[1][target_bin_index]
+
+    # Find the maximum data value that is strictly less than this edge
+    # This mimics: max( intMap[ intMap < h$breaks[index] ] )
+    values_below_edge = filter(x -> x < target_bin_upper_edge, int_values)
+    actual_threshold = isempty(values_below_edge) ? low : maximum(values_below_edge)
+
+    return (low, actual_threshold)
+end
+
+function set_pixel_depth(img, bounds, depth)
+    min_val, max_val = bounds
+    bins = depth - 1
+    
+    if min_val >= max_val
+        return zeros(UInt8, size(img))
+    end
+    
+    # Create intensity bins
+    range_vals = range(min_val, stop=max_val, length=depth)[2:depth]
+    
+    # Assign each pixel to a bin
+    result = similar(img, UInt8) # Use similar to create an array of the same type and size
+    for i in eachindex(img)
+        if img[i] <= min_val
+            result[i] = 0
+        else
+            bin_idx = findfirst(x -> img[i] <= x, range_vals)
+            result[i] = bin_idx === nothing ? bins : bin_idx - 1
+        end
+    end
+    
+    return result
+end
+
+"""
+    norm_slices_hist(slices, bins; prob=0.98)
+
+Normalizes a set of image slices based on a shared histogram range.
+"""
+function norm_slices_hist(slices, bins; prob=0.98)
+    n_files, n_masses = size(slices)
+    norm_img = similar(slices)
+    mass_bounds = []  # This will store bounds for EACH mass
+    
+    # Calculate bounds for each mass across all files
+    for mass_idx in 1:n_masses
+        # Get all slices for this specific mass across all files
+        mass_slices = [slices[i, mass_idx] for i in 1:n_files]
+        all_vals = reduce(vcat, [vec(s) for s in mass_slices])
+        
+        # Calculate global bounds for this specific mass
+        mass_global_bounds = get_outlier_thres(all_vals, prob)
+        push!(mass_bounds, mass_global_bounds)
+        
+        # Normalize each file's slice for this mass using its specific bounds
+        for file_idx in 1:n_files
+            norm_img[file_idx, mass_idx] = set_pixel_depth(slices[file_idx, mass_idx], mass_global_bounds, bins)
+        end
+    end
+    
+    return (norm_img=norm_img, bounds=mass_bounds)  # bounds is now a vector, one per mass
+end
+
+function median_filter(img)
+    # 3x3 median filter implementation
+    return mapwindow(median, img, (3, 3))
+end
+
+# ============================================================================
+#
+#
+# Analysis and Visualization
+#
+# ============================================================================
+
+"""
+    display_statistics(slices)
+
+Calculates and prints key statistics for each slice.
+"""
+function display_statistics(slices, names, masses)
+    if isempty(slices)
+        @warn "Cannot display statistics for empty slice list."
+        return nothing
+    end
+
+    n_files, n_masses = size(slices)
+    #stats_to_calc = Dict("Mean" => mean, "Max" => maximum, "Min" => minimum, "Sum" => sum, "Std" => std)
+    stats_to_calc = Dict("Mean" => mean)
+    all_dfs = Dict{String, DataFrame}()
+
+    for (stat_name, stat_func) in stats_to_calc
+        # Create a matrix to hold the statistic for each slice
+        stat_matrix = zeros(Float64, n_files, n_masses)
+        for i in 1:n_files, j in 1:n_masses
+            flat_slice = vec(slices[i, j])
+            if !isempty(flat_slice)
+                stat_matrix[i, j] = stat_func(flat_slice)
+            end
+        end
+
+        # Create the DataFrame with masses as column headers
+        df = DataFrame(stat_matrix, Symbol.(masses))
+        # Insert the file names as the first column
+        insertcols!(df, 1, :Data => names)
+        mz_row = ["m/z"; masses...]
+        push!(df, mz_row)
+
+        println("\n--- Statistics: $(stat_name) ---")
+        println(df)
+        all_dfs[stat_name] = df
+    end
+    return all_dfs
+end
+
+function plot_slices(slices, names, masses, output_dir; stage_name, bins=256, dpi=150, global_bounds=nothing)
+    n_files, n_masses = size(slices)
+    
+    mkpath(output_dir)
+
+    # If global_bounds is provided, it should now be a VECTOR of bounds (one per mass)
+    # If not provided, calculate per-mass bounds
+    if global_bounds === nothing
+        global_bounds = []
+        for mass_idx in 1:n_masses
+            mass_slices = [slices[i, mass_idx] for i in 1:n_files]
+            all_vals = reduce(vcat, [vec(s) for s in mass_slices])
+            filter!(isfinite, all_vals)
+            mass_bounds = isempty(all_vals) ? (0.0, 1.0) : extrema(all_vals)
+            push!(global_bounds, mass_bounds)
+        end
+    end
+
+    # With Makie, we define a Figure and a layout.
+    fig = Figure(size = (400 * n_masses, 330 * n_files))  # Adjusted size calculation
+
+    # Add a title for the entire figure.
+    Label(fig[0, 1:(2*n_masses)], "Image Slices - $(stage_name)", fontsize=24, font=:bold, tellwidth=false, padding=(0,0,10,0))
+
+    for file_idx in 1:n_files
+        for mass_idx in 1:n_masses
+            img = slices[file_idx, mass_idx]
+            mass = masses[mass_idx]
+            name = names[file_idx]
+            mass_global_bounds = global_bounds[mass_idx]  # Get bounds for THIS specific mass
+
+            # --- Calculate tick properties for THIS mass ---
+            min_val, max_val = mass_global_bounds
+            levels = range(min_val, stop=max_val, length=bins + 1)
+            level_range = levels[end] - levels[1]
+            
+            if level_range == 0
+                levels = range(min_val - 0.1, stop=max_val + 0.1, length=bins + 1)
+                level_range = 0.2
+            end
+            
+            exponent = level_range > 0 ? floor(log10(level_range)) / 3 : 0
+            scale = 10^(3 * exponent)
+            scaled_levels = levels ./ scale
+            
+            format_num = level_range > 0 ? floor(log10(level_range)) % 3 : 0
+            labels = if format_num == 0
+                [@sprintf("%3.2f", lvl) for lvl in scaled_levels]
+            elseif format_num == 1
+                [@sprintf("%3.2f", lvl) for lvl in scaled_levels]
+            else
+                [@sprintf("%3.2f", lvl) for lvl in scaled_levels]
+            end
+            
+            divisors = 2:7
+            remainders = (bins - 1) .% divisors
+            best_divisor = divisors[findlast(x -> x == minimum(remainders), remainders)]
+            tick_indices = round.(Int, range(1, stop=bins + 1, length=best_divisor + 1))
+
+            if !(1 in tick_indices)
+                pushfirst!(tick_indices, 1)
+            end
+            if !((bins + 1) in tick_indices)
+                push!(tick_indices, bins + 1)
+            end
+            unique!(sort!(tick_indices))
+            
+            tick_positions = levels[tick_indices]
+            tick_labels = labels[tick_indices]
+            # --- End of mass-specific tick calculation ---
+
+            # Create an Axis for the heatmap
+            ax = CairoMakie.Axis(fig[file_idx, 2*mass_idx-1],
+                aspect=DataAspect(),
+                title=@sprintf("%s\nm/z: %.2f", basename(name), mass),
+                titlesize=14
+            )
+            hidedecorations!(ax)
+
+            # Use mass-specific bounds for colorrange
+            hm = heatmap!(ax, transpose(img),
+                colormap=cgrad(ColorSchemes.viridis, bins),
+                colorrange=mass_global_bounds  # ← This is mass-specific
+            )
+
+            # Add a colorbar with mass-specific scale
+            cb = Colorbar(fig[file_idx, 2*mass_idx], hm,
+                label=(scale == 1 ? "" : "×10^$(round(Int, 3 * exponent))"),
+                labelpadding=2,
+                labelsize=12,
+                ticks=(tick_positions, tick_labels),
+                ticklabelsize=10
+            )
+            colsize!(fig.layout, 2*mass_idx, 30)
+        end
+    end
+
+    colgap!(fig.layout, 5)
+    rowgap!(fig.layout, 10)
+
+    # Save in multiple formats.
+    formats = ["png", "pdf"]
+    
+    for fmt in formats
+        filename = "$(stage_name)_overview.$(fmt)"
+        save_path = joinpath(output_dir, filename)
+        save(save_path, fig, px_per_unit = dpi / 96.0)
+        @info "Saved $(fmt) overview plot to $save_path"
+    end
+    return fig
+end
+
+# =============================================================================
+#
+#
+# Main Workflow
+#
+# ============================================================================
+
+function load_config(path)
+    config = Dict{String, Any}()
+    if !isfile(path)
+        @warn "Config file not found at $(path). Using default parameters."
+        # Define defaults here in case the file is missing
+        config["masses"] = [100]
+        config["tolerance"] = 0.1
+        config["color_depth"] = 255
+        config["outlier_prob"] = 0.98
+        return config
+    end
+
+    for line in eachline(path)
+        line = strip(line)
+        if isempty(line) || startswith(line, "#")
+            continue # Skip empty lines and comments
+        end
+        parts = split(line, '=', limit=2)
+        if length(parts) != 2
+            @warn "Skipping malformed line in config: $(line)"
+            continue
+        end
+        key = strip(parts[1])
+        value_str = strip(parts[2])
+
+        try
+            if key == "masses"
+                values = [parse(Float64, s) for s in split(value_str, ',')]
+                config[key] = values
+            elseif key == "tolerance"
+                config[key] = parse(Float64, value_str)
+            elseif key == "color_depth"
+                config[key] = parse(Int, value_str)
+            elseif key == "outlier_prob"
+                config[key] = parse(Float64, value_str)
+            end
+        catch e
+            @error "Could not parse value for key '$(key)': $(value_str)"
+        end
+    end
+    return config
+end

src/imzML_old.jl

@@ -0,0 +1,246 @@
+# include( "Common.jl" );
+# *******************************************************************
+# Load Spectra and return a matrix
+#   fileName: Full name path
+# *******************************************************************
+"""
+    LoadImzml( fileName )
+
+Load an imzML file as a matrix. Each column stores x-pixel position,
+ y-pixel position, x-axis data and y-axis data.
+
+# Arguments
+* `fileName`: Full path name of the imzML file
+
+# Examples
+```julia
+# Load DESI MSI Carcinoma image data
+spectra = LoadImzml( "80TopL, 50TopR, 70BottomL, 60BottomR-centroid.imzML" )
+size( spectra )
+(4, 18632)
+```
+"""
+function LoadImzml( fileName )
+
+  # Open file handles
+   if !isfile(fileName)
+        error("provided path is not a file")
+    end
+
+    if endswith(fileName, ".imzML")
+        stream = open(fileName)
+        hIbd = open(replace(fileName, ".imzML" => ".ibd"))
+    else
+        stream = open(fileName * ".imzML")
+        hIbd = open(fileName * ".ibd")
+    end
+
+  # Get axes types and image dimensions
+  axis   = AxesConfigImg( stream )
+  imgDim = GetImgDimensions( stream )
+  format = [ axis[1].Format, axis[2].Format ]
+
+  # Locate spectrum attributes
+  start  = position( stream )  
+  attr   = GetSpectrumAttributes( stream, hIbd ) 
+
+  # Load spectra
+  seek( stream, start )
+  spectra = LoadImgData( stream, hIbd, attr, imgDim[3], format )
+
+  close( stream )
+  close( hIbd )
+  return spectra
+
+end
+
+
+# *******************************************************************
+# Get axes value type
+# *******************************************************************
+function AxesConfigImg( stream )
+
+  # Locate which axis is defined at first
+  tag   = FindTag( stream, r"^\s*<(referenceableParamGroup )" )
+  value = GetAttribute( tag.captures[1], "intensityArray" )
+  order = 1 + ( value !== nothing )
+
+  # Read first axis configuration
+  axis  = Array{ SpecDim, 1 }( undef, 2 )
+  axis[ order ] = ConfigureSpecDim( stream )
+
+  # Read second axis configuration
+  FindTag( stream, r"^\s*<(referenceableParamGroup )" )
+  axis[ xor( order,3 ) ] = ConfigureSpecDim( stream )
+  return axis
+
+end
+
+
+# *******************************************************************
+# Get vector's storage options and image dimensions
+# *******************************************************************
+function GetImgDimensions( stream )
+
+  # Looks for "scanSettings" tag
+  # FindImgTag( stream, "scanSettings" )
+  FindTag( stream, r"^\s*<(scanSettings )" )
+
+  # Initial values for dimension retrieve
+  n         = 2
+  dim       = [ 0, 0, 0 ]
+  currLine  = ""
+  matchInfo = RegexMatch
+  
+  while( n > 0 )
+  
+    # Next field
+    currLine  = readline( stream )
+    matchInfo = match( r"^\s*<(cvParam)", currLine )
+    index     = length( matchInfo.captures[1] )
+    matchInfo = GetAttribute( currLine[index:end], "accession" )
+    
+    # Get axis identity
+    if( matchInfo.captures[1] == "IMS:1000042"        # max X
+    ||  matchInfo.captures[1] == "IMS:1000043" )      # max Y
+    
+      # Read dimension's pixels 
+      axis      = matchInfo.captures[1][end] - '1'
+      index    += matchInfo.offsets[1] + length( matchInfo.captures[1] )
+      matchInfo = GetAttribute( currLine[index:end], "value" )
+      dim[axis] = parse( Int32, matchInfo.captures[1] )
+      n        -= 1
+    end
+  end
+
+  # Load stored spectra counter
+  matchInfo = FindTag( stream, r"^\s*<spectrumList(.+)" )
+  matchInfo = GetAttribute( matchInfo.captures[1], "count" )
+  dim[ 3 ]  = parse( Int32, matchInfo.captures[1] )
+  
+  return dim
+end
+
+
+# *******************************************************************
+# Length of "spectrum" tag without attribute values
+# *******************************************************************
+function  GetSpectrumTag( stream )
+
+  offset = position( stream )
+  tag    = FindTag( stream, r"^\s*<spectrum (.+)" )
+  first  = 1
+
+  while true
+    
+    value = match( r"[^=]+=\"([^\"]+)\"", tag.captures[1][first:end] )
+    if value === nothing
+      break
+    end
+    first  += value.offsets[1] + length( value.captures[1] )
+    offset += length( value.captures[1] )
+  end
+
+  return offset
+end
+
+# seek( stream, 7141 )
+
+# *******************************************************************
+# Locate spectrum attribute's 
+#   [1]  1: x position stored first   2: y position stored first
+#   [2]  Second stored position
+#   [3]  1: mzArray stored first      2: intensityArray stored first
+#   [4]  Second stored axis
+#   [5]  skip chars to find first dimension
+#   [6]  skip chars to find second dimension
+#   [7]  skip chars to find vector length
+#   [8]  skip chars to find next spectrum  
+# *******************************************************************
+function GetSpectrumAttributes( stream, hIbd )
+
+  # Reserve memory for vector configuration
+  skip   = Vector{ UInt32 }( undef, 8 )
+  offset = GetSpectrumTag( stream )
+
+  # Get axis order and position of pixel coordinate values
+  tag     = FindTag( stream, r" accession=\"IMS:100005(\d)\"(.+)" )
+  skip[1] = tag.captures[1][1] + 1 - '0'
+  skip[2] = xor( skip[1], 3 )
+
+  value   = match( r" value=\"(\d+)\".+", tag.captures[2] )
+  skip[5] = position( stream ) - offset - length( value.match ) - 2
+  value   = match( r" value=\"\d+\".+", readline( stream ) )
+  skip[6] = value.offset
+
+  # Get axis order 
+  offset  = position( stream )
+  tag     = FindTag( stream, r"^\s*<referenceableParamGroupRef(.+)" )
+  value   = GetAttribute( tag.captures[1], "ref" )
+  skip[3] = ( value.captures[1] == "intensityArray" ) + 3
+  skip[4] = xor( skip[3], 7 )
+
+  k = 2
+  while k != 0
+
+    tag = FindTag( stream, r" accession=\"IMS:100010(\d)\"(.+)" )
+
+    # Set IBD offset for first spectrum
+    if tag.captures[1][1] == '2'
+      value = GetAttribute( tag.match, "value" )
+      seek( hIbd, parse( Int32, value.captures[1] ) )
+      k -= 1
+      continue
+
+    elseif tag.captures[1][1] == '3'
+      skip[7]  = position( stream ) - offset - length( tag.match )
+      offset   = position( stream )
+      k -= 1
+    end
+  end
+
+  # Compute characters to skip 
+  FindTag( stream, r"^\s*</spectrum" )
+  skip[8] = position( stream ) - offset
+  return skip
+
+end
+
+
+# *******************************************************************
+# Read spectral data
+# *******************************************************************
+function LoadImgData( stream, hIbd, attr, imgDim, format )
+
+  # Reserve spectra memory & update IBD seek offset
+  spectra = Array{Any}( undef, ( 4, imgDim ) )
+  contador = 0;
+
+  for k in 1:imgDim
+
+    # Load image coordinates
+    skip( stream, attr[5] )
+    value                = FindTag( stream, r"value=\"(\d+)\"" )
+    spectra[ attr[1],k ] = parse( Int32, value.captures[1] )
+
+    skip( stream, attr[6] )
+    value                = FindTag( stream, r"value=\"(\d+)\"" )
+    spectra[ attr[2],k ] = parse( Int32, value.captures[1] )
+
+    # Get vector length
+    skip( stream, attr[7] )
+    value   = FindTag( stream, r"value=\"(\d+)\"" )
+    nPoints = parse( Int32, value.captures[1] )
+
+    # Reserve memory and read spectrum
+    spectra[ 3,k ] = Array{ format[1] }( undef, nPoints )
+    spectra[ 4,k ] = Array{ format[2] }( undef, nPoints )
+    read!( hIbd, spectra[ attr[3],k ] )
+    read!( hIbd, spectra[ attr[4],k ] )
+    skip( stream, attr[8] )
+
+  end
+
+  return spectra
+
+end

src/mzML.jl

@@ -0,0 +1,191 @@
+# /home/pixel/Documents/Cinvestav_2025/JuliaMSI/src/mzML_lazy.jl
+
+using Base64, Libz
+include("ParserHelpers.jl")
+
+# ============================================================================
+# Data Structures for Lazy Loading
+# ============================================================================ 
+
+struct SpectrumAsset
+    format::Type
+    is_compressed::Bool
+    offset::Int64
+    encoded_length::Int32
+    axis_type::Symbol # :mz or :intensity
+end
+
+struct MzMLSpectrumMetadata
+    id::String
+    mz_asset::SpectrumAsset
+    intensity_asset::SpectrumAsset
+end
+
+mutable struct MzMLData
+    file_handle::IO
+    spectra_metadata::Vector{MzMLSpectrumMetadata}
+
+    function MzMLData(file_handle::IO, spectra_metadata::Vector{MzMLSpectrumMetadata})
+        obj = new(file_handle, spectra_metadata)
+        finalizer(obj) do o
+            if isopen(o.file_handle)
+                close(o.file_handle)
+            end
+        end
+        return obj
+    end
+end
+
+# ============================================================================ 
+# mzML Parser Implementation
+# ============================================================================ 
+
+function get_spectrum_asset_metadata(stream)
+    start_pos = position(stream)
+    
+    bda_tag = find_tag(stream, r"<binaryDataArray\s+encodedLength=\"(\d+)\"")
+    if bda_tag === nothing; error("Cannot find binaryDataArray"); end
+    encoded_length = parse(Int32, bda_tag.captures[1])
+
+    params = CVParams(Float64, false, :mz)
+
+    while true
+        line = readline(stream)
+        if occursin("</binaryDataArray>", line); break; end
+        if occursin("<cvParam", line)
+            accession = get_attribute(line, "accession")
+            if accession === nothing; continue; end
+            update_cv_params!(params, accession.captures[1])
+        end
+    end
+
+    seek(stream, start_pos)
+    binary_tag = find_tag(stream, r"<binary>")
+    if binary_tag === nothing; error("Cannot find binary tag"); end
+    binary_offset = position(stream)
+
+    find_tag(stream, r"</binaryDataArray>") # Position for next call
+
+    return SpectrumAsset(params.format, params.is_compressed, binary_offset, encoded_length, params.axis_type)
+end
+
+function parse_spectrum_metadata(stream, offset::Int64)
+    seek(stream, offset)
+    
+    id_match = find_tag(stream, r"<spectrum\s+index=\"\d+\"\s+id=\"([^"]+)\"")
+    id = id_match === nothing ? "" : id_match.captures[1]
+
+    asset1 = get_spectrum_asset_metadata(stream)
+    asset2 = get_spectrum_asset_metadata(stream)
+
+    mz_asset = asset1.axis_type == :mz ? asset1 : asset2
+    int_asset = asset1.axis_type == :intensity ? asset1 : asset2
+
+    return MzMLSpectrumMetadata(id, mz_asset, int_asset)
+end
+
+function parse_offset_list(stream)
+    offsets = Int64[]
+    while true
+        line = readline(stream)
+        m = match(r"<offset\s+idRef=\"\S+\">(\d+)</offset>", line)
+        if m !== nothing
+            push!(offsets, parse(Int64, m.captures[1]))
+        elseif occursin("</index>", line) || occursin("</indexedmzML>", line)
+            break
+        end
+    end
+    return offsets
+end
+
+function load_mzml_lazy(file_path::String)
+    stream = open(file_path)
+
+    try
+        seekend(stream)
+        skip(stream, -4096) # Read last 4KB
+        footer = read(stream, String)
+        
+        index_offset_match = match(r"<indexListOffset>(\d+)</indexListOffset>", footer)
+        if index_offset_match === nothing
+            error("Could not find <indexListOffset>. File may not be an indexed mzML.")
+        end
+        
+        index_offset = parse(Int64, index_offset_match.captures[1])
+        seek(stream, index_offset)
+        
+        if find_tag(stream, r"<index\s+name=\"spectrum\"") === nothing
+            error("Could not find spectrum index.")
+        end
+
+        spectrum_offsets = parse_offset_list(stream)
+        if isempty(spectrum_offsets)
+            error("No spectrum offsets found.")
+        end
+
+        spectra_metadata = [parse_spectrum_metadata(stream, offset) for offset in spectrum_offsets]
+        
+        # We return the stream inside MzMLData, so we don't close it here.
+        # The finalizer on MzMLData will handle it.
+        return MzMLData(stream, spectra_metadata)
+
+    catch e
+        close(stream) # Close file on error
+        rethrow(e)
+    end
+end
+
+# ============================================================================ 
+# Data Access Functions
+# ============================================================================ 
+
+function read_binary_vector(io::IO, asset::SpectrumAsset)
+    seek(io, asset.offset)
+    
+    # Read the raw base64 data. The regex is slow. Reading bytes is better.
+    # The old code did `Vector{UInt8}(base64Vec.captures[1])`. This implies ASCII.
+    raw_b64 = read(io, asset.encoded_length)
+    
+    decoded_bytes = base64decode(String(raw_b64))
+    
+    bytes_io = IOBuffer(asset.is_compressed ? Libz.inflate(decoded_bytes) : decoded_bytes)
+    
+    n_elements = Int(bytes_io.size / sizeof(asset.format))
+    out_array = Array{asset.format}(undef, n_elements)
+    read!(bytes_io, out_array)
+    
+    return out_array
+end
+
+function get_spectrum(data::MzMLData, index::Int)
+    meta = data.spectra_metadata[index]
+    
+    mz_array = read_binary_vector(data.file_handle, meta.mz_asset)
+    intensity_array = read_binary_vector(data.file_handle, meta.intensity_asset)
+    
+    return (id=meta.id, mz=mz_array, intensity=intensity_array)
+end
+
+"""
+    LoadMzml(fileName::String)
+
+Eagerly loads all spectra from a .mzML file.
+Provided for backward compatibility.
+"""
+function LoadMzml(fileName::String)
+    mzml_data = load_mzml_lazy(fileName)
+    
+    num_spectra = length(mzml_data.spectra_metadata)
+    spectra_matrix = Array{Any}(undef, (2, num_spectra))
+    
+    for i in 1:num_spectra
+        spectrum = get_spectrum(mzml_data, i)
+        spectra_matrix[1, i] = spectrum.mz
+        spectra_matrix[2, i] = spectrum.intensity
+    end
+    
+    # Eager load, so we can close the handle now.
+    close(mzml_data.file_handle)
+    
+    return spectra_matrix
+end

src/mzML_old.jl

@@ -0,0 +1,350 @@
+# *******************************************************************
+# Load scans from mzML file
+#   fileName: Full name path
+# *******************************************************************
+"""
+    `LoadMzml( fileName )`
+
+Load the mzML file as a Matrix. Each column stores a single scan,
+ x-axis and y-axis scan data are stored in the first and second
+ row respectively.
+
+# Arguments
+* `fileName`: Full path name of the mzML file
+
+# Examples
+```julia
+spectra = LoadMzml( "T9_A1.mzML" )
+size( spectra )
+(2, 10)
+using Plots
+plot( spectra[1,4], spectra[2,4] )  
+```
+
+The following image will appear in an auxiliary window
+
+![](./Plots.bmp)
+"""
+function LoadMzml( fileName )
+
+  # Open file hFile = open( joinpath( data_dir, "Col_1.mzML" ), "r+" )
+  stream    = open( fileName, "r+" )
+	intensity = []
+
+  # Load footer
+  seekend( stream )
+  skip( stream, -160 )
+  raw    = read( stream, 160 )
+  footer = String( raw )
+
+  # Search "indexListOffset" tag
+  index = match( r"<indexListOffset>(?<Ofst>\d+)</indexListOffset>", footer )
+
+  if index !== nothing
+
+    # Jump to index table
+    seek( stream, parse( Int64, index.captures[1] ) )
+
+		# Find first "index" tag
+		index = FindTag( stream, r"<index +name=\"spectrum" )
+
+		# Load offset list and load spectra
+		if index !== nothing
+			intensity = LoadSpectra( stream, ParseOffsetList( stream ) )
+		end
+
+	end
+
+  # Close file handle
+  close( stream )
+
+  # Return spectral matrix
+  return intensity
+
+end
+
+
+# *******************************************************************
+# Parse indexList offsets
+#   stream: Source file stream
+# *******************************************************************
+function ParseOffsetList( stream )
+
+  # List of 8 elements
+  entry  = zeros( Int64, 8 )
+  used   = 0
+  unused = 8
+
+  while true
+
+    # Get offset string
+    raw   = readline( stream )
+    index = match( r">(?<Ofst>\d+)<", String( raw ) )
+
+    # Exit condition
+    if index === nothing
+      break
+    end
+
+    # Save index
+    used   += 1
+    unused -= 1
+    entry[used] = parse( Int64, index.captures[1] )
+
+    # Grow list when needed
+    if unused == 0
+      append!( entry, zeros( Int64, 8 ) )
+      unused = 8
+    end
+
+  end
+
+  return entry[1:used]
+
+end
+
+
+# *******************************************************************
+# Read vector data
+#   stream: Source file stream
+# *******************************************************************
+function ReadVector( stream, axis )
+
+	# Read base64 vector
+	base64Vec = FindTag( stream, r"<binary>([^<]+)<" )
+  base64Vec = Vector{ UInt8 }( base64Vec.captures[1] )
+
+  # Fill a IO stream with decoded data
+	if axis.Packed == 1
+		io = IOBuffer( Libz.inflate( Decode64( base64Vec ) ) )
+	else
+		io = IOBuffer( Decode64( base64Vec ) )
+	end
+
+  # Convert IO stream content as dataType vector
+  nElem = Int32( io.size / sizeof( axis.Format ) )
+  out   = Array{ axis.Format }( undef, nElem )
+  read!( io, out )
+
+	return out
+
+end
+
+
+# *******************************************************************
+# Get the decoding options for the first vector in spectrum
+#    stream: Valid file stream
+#   lastTag: Last tag to skip i.e. "binaryDataArray"
+# *******************************************************************
+function AxisConfig( stream )
+
+  offset = position( stream )
+  skip   = 0
+  while true
+
+    # Load line
+    raw   = readline( stream )
+    tag   = match( r"^\s*<([^\s]+)", raw )
+    start = length( tag.match ) + 1
+
+    # Is a "cvParam" tag?
+    if tag.captures[1] == "cvParam"
+
+      # Get accesion field value
+      attribute = GetAttribute( raw[ start:end ], "accession" )
+
+      # Subtract file name length from skip field
+      if attribute.captures[1] != "MS:1000796"
+        continue
+      else
+        start += attribute.offset + length( attribute.match )
+        value  = GetAttribute( raw[ start:end ], "value" )
+        skip  += length( value.captures[1] )
+        continue
+      end
+
+    # Exit condition
+    elseif tag.captures[1] == "binaryDataArray"
+      value  = GetAttribute( raw[ start:end ], "encodedLength" )
+      offset = position( stream ) - offset
+      front  = ConfigureSpecDim( stream )
+      front.Skip += offset - length( value.captures[1] ) - skip
+      return ( front )
+
+    end
+  end
+end
+
+
+# *******************************************************************
+# Load Spectra and return a matrix
+#   stream: Source file stream
+# *******************************************************************
+function LoadSpectra( stream, offset )
+
+  # Move file pointer towards first spectrum tag
+  seek( stream, offset[1] )
+
+  # Decodes axis information and skip counter
+  firstAxis  = AxisConfig( stream )
+  secondAxis = AxisConfig( stream )
+
+  # Allocates memory for spectra storage
+	specCount = length( offset )
+	spectrum  = Array{ Any }( undef, ( 2, specCount ) )
+	index     = 1
+  order     = 1 + ( firstAxis.Axis == 2 )
+
+	for i in 1:specCount
+
+		# Load first vector
+    seek( stream, offset[ i ] )
+    skip( stream, firstAxis.Skip )
+		spectrum[ order, i ] = ReadVector( stream, firstAxis )
+
+    # Load second vector
+    skip( stream, secondAxis.Skip )
+		spectrum[ xor( order,3 ), i ] = ReadVector( stream, secondAxis )
+	end
+
+  return( spectrum )
+
+end
+
+
+# *******************************************************************
+# Decoding base64 blocks
+# *******************************************************************
+
+# Decoding exchange table
+base64dec = Vector{UInt8}( [
+  62, # (01) + -> 43d, 0x2B
+   0, # (02) ,
+   0, # (03) -
+   0, # (04) .
+  63, # (05) /
+  52, # (06) 0
+  53, # (07) 1
+  54, # (08) 2
+  55, # (09) 3
+  56, # (10) 4
+  57, # (11) 5
+  58, # (12) 6
+  59, # (13) 7
+  60, # (14) 8
+  61, # (15) 9
+   0, # (16) :
+   0, # (17) ;
+   0, # (18) <
+   0, # (19) =
+   0, # (20) >
+   0, # (21) ?
+   0, # (22) @
+   0, # (23) A
+   1, # (24) B
+   2, # (25) C
+   3, # (26) D
+   4, # (27) E
+   5, # (28) F
+   6, # (29) G
+   7, # (40) H
+   8, # (41) I
+   9, # (42) J
+  10, # (43) K
+  11, # (44) L
+  12, # (45) M
+  13, # (46) N
+  14, # (47) O
+  15, # (48) P
+  16, # (49) Q
+  17, # (40) R
+  18, # (41) S
+  19, # (42) T
+  20, # (43) U
+  21, # (44) V
+  22, # (45) W
+  23, # (46) X
+  24, # (47) Y
+  25, # (48) Z
+   0, # (49) [
+   0, # (50) \
+   0, # (51) ]
+   0, # (52) ^
+   0, # (53) _
+   0, # (54) `
+  26, # (55) a
+  27, # (56) b
+  28, # (57) c
+  29, # (58) d
+  30, # (59) e
+  31, # (60) f
+  32, # (61) g
+  33, # (62) h
+  34, # (63) i
+  35, # (64) j
+  36, # (65) k
+  37, # (66) l
+  38, # (67) m
+  39, # (68) n
+  40, # (69) o
+  41, # (70) p
+  42, # (71) q
+  43, # (72) r
+  44, # (73) s
+  45, # (74) t
+  46, # (75) u
+  47, # (76) v
+  48, # (77) w
+  49, # (78) x
+  50, # (79) y
+  51, # (80) z
+] )
+
+# base64 decoding
+function DecodeTriplet( data, base64 )
+
+  # Character to binary equivalent
+  data    = data .- 0x2A
+  data[1] = base64[ data[1] ]
+  data[2] = base64[ data[2] ]
+  data[3] = base64[ data[3] ]
+  data[4] = base64[ data[4] ]
+
+  # Get Triplet
+  decode    = Array{UInt8}( undef,3 )
+  decode[1] = data[1] << 2 + data[2] >> 4
+  decode[2] = data[2] << 4 + data[3] >> 2
+  decode[3] = data[3] << 6 + data[4]
+  return decode
+
+end
+
+
+function Decode64( data )
+
+  # Reserve uninitialized memory
+  block  = divrem( sizeof( data ), 4 )
+  decode = Array{UInt8}( undef, 3*block[1] )
+
+  # Recover binary from Chars
+  index = 1
+  first = 1
+  last  = 3
+
+  for i in 1:block[1]
+    decode[first:last] = DecodeTriplet( data[index:index+3], base64dec )
+    index += 4
+    first += 3
+    last  += 3
+  end
+
+  # Adjust block size
+  if data[end-1] == 0x3D
+    decode = decode[1:end-2]
+  elseif data[end] == 0x3D
+    decode = decode[1:end-1]
+  end
+
+  return decode
+
+end

start_MSI_GUI.jl

@@ -1,7 +1,11 @@
 using Pkg
 sTime=time()
 Pkg.activate(".")
-Pkg.instantiate()
+# Only instantiate in development mode
+if get(ENV, "GENIE_ENV", "dev") != "prod"
+    @info "Development environment detected. Instantiating packages..."
+    Pkg.instantiate()
+end
 Pkg.gc()
 
 using Genie