From b1d7838a6002c3f9293c30a54a039306ec6023c8 Mon Sep 17 00:00:00 2001 From: Pixelguy14 Date: Thu, 25 Sep 2025 10:42:33 -0600 Subject: [PATCH] Implemented imzml and mzml lazy loading, imzml to load spectra matrix, working on mzml to imzml --- .gitignore | 1 + app.jl | 33 +- config/env/global.jl | 3 +- julia_imzML_visual.jl | 2 +- "ongoing_tasks.txt\n" | 27 -- scripts/build.jl | 23 ++ scripts/precompile.jl | 45 +++ src/Bitmap.jl | 399 ++++++++++++++++++++++++ src/Common.jl | 94 ++++++ src/DataAccess.jl | 82 +++++ src/MzmlConverter.jl | 270 ++++++++++++++++ src/ParserHelpers.jl | 122 ++++++++ src/imzML.jl | 699 ++++++++++++++++++++++++++++++++++++++++++ src/imzML_old.jl | 246 +++++++++++++++ src/mzML.jl | 191 ++++++++++++ src/mzML_old.jl | 350 +++++++++++++++++++++ start_MSI_GUI.jl | 6 +- 17 files changed, 2547 insertions(+), 46 deletions(-) delete mode 100644 "ongoing_tasks.txt\n" create mode 100644 scripts/build.jl create mode 100644 scripts/precompile.jl create mode 100644 src/Bitmap.jl create mode 100644 src/Common.jl create mode 100644 src/DataAccess.jl create mode 100644 src/MzmlConverter.jl create mode 100644 src/ParserHelpers.jl create mode 100644 src/imzML.jl create mode 100644 src/imzML_old.jl create mode 100644 src/mzML.jl create mode 100644 src/mzML_old.jl diff --git a/.gitignore b/.gitignore index 6832527..1905238 100644 --- a/.gitignore +++ b/.gitignore @@ -4,3 +4,4 @@ public/css/imgOver.png !public/css/autogenerated.css !public/css/LABI_logo.png log/* +R_original_scripts/ diff --git a/app.jl b/app.jl index aaa9767..4bf8124 100644 --- a/app.jl +++ b/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 ySpectraMz=Float64[] + @out xSpectraMz = Vector{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 diff --git a/config/env/global.jl b/config/env/global.jl index 0b24049..94f0b9d 100644 --- a/config/env/global.jl +++ b/config/env/global.jl @@ -1 +1,2 @@ -ENV["GENIE_ENV"] = "prod" +ENV["GENIE_ENV"] = "dev" +#ENV["GENIE_ENV"] = "prod" diff --git a/julia_imzML_visual.jl b/julia_imzML_visual.jl index 1316186..604fcb1 100644 --- a/julia_imzML_visual.jl +++ b/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 diff --git "a/ongoing_tasks.txt\n" "b/ongoing_tasks.txt\n" deleted file mode 100644 index 94a226e..0000000 --- "a/ongoing_tasks.txt\n" +++ /dev/null @@ -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 diff --git a/scripts/build.jl b/scripts/build.jl new file mode 100644 index 0000000..86eafa9 --- /dev/null +++ b/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" diff --git a/scripts/precompile.jl b/scripts/precompile.jl new file mode 100644 index 0000000..48b33da --- /dev/null +++ b/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." diff --git a/src/Bitmap.jl b/src/Bitmap.jl new file mode 100644 index 0000000..a7edbf1 --- /dev/null +++ b/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 ] + diff --git a/src/Common.jl b/src/Common.jl new file mode 100644 index 0000000..aa26b94 --- /dev/null +++ b/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 + + diff --git a/src/DataAccess.jl b/src/DataAccess.jl new file mode 100644 index 0000000..5977248 --- /dev/null +++ b/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 diff --git a/src/MzmlConverter.jl b/src/MzmlConverter.jl new file mode 100644 index 0000000..16a39ba --- /dev/null +++ b/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" 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>", 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" 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 + # - , , + # - A with a for each pixel. + # - Each 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 diff --git a/src/ParserHelpers.jl b/src/ParserHelpers.jl new file mode 100644 index 0000000..af53599 --- /dev/null +++ b/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 + + diff --git a/src/imzML.jl b/src/imzML.jl new file mode 100644 index 0000000..2c87316 --- /dev/null +++ b/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("` tag, ignoring attribute values. +""" +function get_spectrum_tag_offset(stream) + offset = position(stream) + tag = find_tag(stream, r"^\s*") + 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 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 diff --git a/src/imzML_old.jl b/src/imzML_old.jl new file mode 100644 index 0000000..dbb1480 --- /dev/null +++ b/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*", line); break; end + if occursin("") + if binary_tag === nothing; error("Cannot find binary tag"); end + binary_offset = position(stream) + + find_tag(stream, r"") # 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"(\d+)", line) + if m !== nothing + push!(offsets, parse(Int64, m.captures[1])) + elseif occursin("", line) || occursin("", 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"(\d+)", footer) + if index_offset_match === nothing + error("Could not find . 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"(?\d+)", footer ) + + if index !== nothing + + # Jump to index table + seek( stream, parse( Int64, index.captures[1] ) ) + + # Find first "index" tag + index = FindTag( stream, r"(?\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"([^<]+)<" ) + 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 diff --git a/start_MSI_GUI.jl b/start_MSI_GUI.jl index 0d69ef8..14036e2 100644 --- a/start_MSI_GUI.jl +++ b/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