-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
+
+
+
+
+
+
+
+ Image visualizer
+
+
+
+
+
-
-
-
-
-
+
+
+
+
+
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
+
+
-
-
-
-
-
+
+
+ TrIQ visualizer
+
+
+
+
+
+
+
+
+
+
-
-
-
- Batch Process Summary
-
+
+
+
+
+
+
+ Loading plot
+
-
- {{ batch_summary }}
-
-
-
-
-
-
-
-
-
-
-
- Dataset Summary
-
-
-
-
-
-
-
-
- {{ row.parameter }}
-
-
- {{ row.value }}
-
+
+
+
+ Mean spectrum plot
+
-
-
-
-
No metadata found in registry for this dataset.
-
You may need to process the file first.
-
+
\ No newline at end of file
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+ Warning
+
+
+
+ {{msg}}
+
+
+
+
+
+
+
+
+
+
+
+ Compare two different images or plots
+
+
+ Transparency: {{ imgTrans }}
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+ Batch Process Summary
+
+
+
+ {{ batch_summary }}
+
+
+
+
+
+
+
+
+
+
+
+ Dataset Summary
+
+
+
+
+
+
+
+
+ {{ row.parameter }}
+
+
+ {{ row.value }}
+
+
+
+
+
+
No metadata found in registry for this dataset.
+
You may need to process the file first.
+
\ No newline at end of file
diff --git a/julia_imzML_visual.jl b/julia_imzML_visual.jl
index aff0afa..803af8d 100644
--- a/julia_imzML_visual.jl
+++ b/julia_imzML_visual.jl
@@ -1,32 +1,67 @@
-# == Search functions ==
-# Functions that recieve a list to update, and the current direction both as string for
-# searching in the directory the position the list is going
+# julia_imzML_visual.jl
+
+"""
+ increment_image(current_image, image_list)
+
+Finds the next image in a list.
+
+# Arguments
+- `current_image`: The current image file name.
+- `image_list`: The list of available image file names.
+
+# Returns
+- The file name of the next image, or the last image if the current one is the last or not found.
+- `nothing` if `image_list` is empty.
+"""
function increment_image(current_image, image_list)
if isempty(image_list)
return nothing
end
- current_index=findfirst(isequal(current_image), image_list)
- if current_index==nothing || current_index==length(image_list) || current_image ===""
- return image_list[length(image_list)] # Return the current image if it's the last one or not found
+ current_index = findfirst(isequal(current_image), image_list)
+ if current_index === nothing || current_index == length(image_list) || current_image == ""
+ return image_list[end] # Return the last image if current is not found or is the last
else
return image_list[current_index + 1] # Move to the next image
end
end
+"""
+ decrement_image(current_image, image_list)
+
+Finds the previous image in a list.
+
+# Arguments
+- `current_image`: The current image file name.
+- `image_list`: The list of available image file names.
+
+# Returns
+- The file name of the previous image, or the first image if the current one is the first or not found.
+- `nothing` if `image_list` is empty.
+"""
function decrement_image(current_image, image_list)
if isempty(image_list)
return nothing
end
- current_index=findfirst(isequal(current_image), image_list)
- if current_index==nothing || current_index==1 || current_image===""
- return image_list[1] # Return the current image if it's the first one or not found
+ current_index = findfirst(isequal(current_image), image_list)
+ if current_index === nothing || current_index == 1 || current_image == ""
+ return image_list[1] # Return the first image if current is not found or is the first
else
return image_list[current_index - 1] # Move to the previous image
end
end
-## Plot Image functions
-# Downsample an image matrix to a maximum dimension while preserving aspect ratio
+"""
+ downsample_image(img_matrix, max_dim::Int)
+
+Downsamples an image matrix to a maximum dimension while preserving aspect ratio.
+
+# Arguments
+- `img_matrix`: The image matrix to downsample.
+- `max_dim`: The maximum dimension (width or height) for the downsampled image.
+
+# Returns
+- The downsampled image matrix.
+"""
function downsample_image(img_matrix, max_dim::Int)
h, w = size(img_matrix)
if h <= max_dim && w <= max_dim
@@ -46,26 +81,37 @@ function downsample_image(img_matrix, max_dim::Int)
return imresize(img_matrix, (new_h, new_w))
end
-# loadImgPlot recieves the local directory of the image as a string,
-# returns the layout and data for the heatmap plotly plot
-# this function loads the image into a plot
+"""
+ loadImgPlot(interfaceImg::String)
+
+Loads an image and creates a Plotly heatmap.
+
+# Arguments
+- `interfaceImg`: The path to the image file relative to the "public" directory.
+
+# Returns
+- `plotdata`: A vector containing the Plotly trace.
+- `plotlayout`: The Plotly layout for the plot.
+- `width`: The width of the loaded image.
+- `height`: The height of the loaded image.
+"""
function loadImgPlot(interfaceImg::String)
# Load the image
- cleaned_img=replace(interfaceImg, r"\?.*" => "")
- cleaned_img=lstrip(cleaned_img, '/')
- var=joinpath("./public", cleaned_img)
- img=load(var)
+ cleaned_img = replace(interfaceImg, r"\?.*" => "")
+ cleaned_img = lstrip(cleaned_img, '/')
+ var = joinpath("./public", cleaned_img)
+ img = load(var)
# Convert to grayscale
- img_gray=Gray.(img)
- img_array=Array(img_gray)
- elevation=Float32.(Array(img_array)) ./ 255.0
+ img_gray = Gray.(img)
+ img_array = Array(img_gray)
+ elevation = Float32.(Array(img_array)) ./ 255.0
# Get the X, Y coordinates of the image
- height, width=size(img_array)
- X=collect(1:width)
- Y=collect(1:height)
+ height, width = size(img_array)
+ X = 1:width
+ Y = 1:height
# Create the layout
- layout=PlotlyBase.Layout(
+ layout = PlotlyBase.Layout(
title=PlotlyBase.attr(
text="",
font=PlotlyBase.attr(
@@ -83,11 +129,11 @@ function loadImgPlot(interfaceImg::String)
visible=false,
range=[-height, 0]
),
- margin=attr(l=0,r=0,t=0,b=0,pad=0)
+ margin=attr(l=0, r=0, t=0, b=0, pad=0)
)
# Create the trace for the image
- trace=PlotlyBase.heatmap(
+ trace = PlotlyBase.heatmap(
z=elevation,
x=X,
y=-Y,
@@ -113,16 +159,29 @@ function loadImgPlot(interfaceImg::String)
)
)
- plotdata=[trace]
- plotlayout=layout
+ plotdata = [trace]
+ plotlayout = layout
return plotdata, plotlayout, width, height
end
-# loadImgPlot recieves the local directory of the image as a string, the local directory o the overlay image
-# and the transparency its required to have. Returns the layout and data for the heatmap plotly plot
-# this function loads the image into a plot
+"""
+ loadImgPlot(interfaceImg::String, overlayImg::String, imgTrans::Float64)
+
+Loads a main image and overlays a second image on top, creating a Plotly heatmap.
+
+# Arguments
+- `interfaceImg`: Path to the main image file.
+- `overlayImg`: Path to the overlay image file.
+- `imgTrans`: Transparency level for the overlay image (0.0 to 1.0).
+
+# Returns
+- `plotdata`: A vector containing the Plotly trace for the main image.
+- `plotlayout`: The Plotly layout, including the overlay image.
+- `width`: The width of the main image.
+- `height`: The height of the main image.
+"""
function loadImgPlot(interfaceImg::String, overlayImg::String, imgTrans::Float64)
- timestamp=string(time_ns())
+ timestamp = string(time_ns())
# Load the main image
cleaned_img = replace(interfaceImg, r"\?.*" => "")
cleaned_img = lstrip(cleaned_img, '/')
@@ -134,8 +193,8 @@ function loadImgPlot(interfaceImg::String, overlayImg::String, imgTrans::Float64
elevation = Float32.(Array(img_array)) ./ 255.0
# Get the X, Y coordinates of the image
height, width = size(img_array)
- X = collect(1:width)
- Y = collect(1:height)
+ X = 1:width
+ Y = 1:height
# Create the layout with overlay image
layoutImg = PlotlyBase.Layout(
@@ -147,40 +206,40 @@ function loadImgPlot(interfaceImg::String, overlayImg::String, imgTrans::Float64
color="black"
)
),
- images = [attr(
- source = "$(overlayImg)?t=$(timestamp)",
- xref = "x",
- yref = "y",
- x = 0,
- y = 0,
- sizex = width,
- sizey = -height,
- sizing = "stretch",
- opacity = imgTrans,
- layer = "above" # Place the overlay image in the foreground
+ images=[attr(
+ source="$(overlayImg)?t=$(timestamp)",
+ xref="x",
+ yref="y",
+ x=0,
+ y=0,
+ sizex=width,
+ sizey=-height,
+ sizing="stretch",
+ opacity=imgTrans,
+ layer="above" # Place the overlay image in the foreground
)],
- xaxis = PlotlyBase.attr(
- visible = false,
- scaleanchor = "y",
- range = [0, width]
+ xaxis=PlotlyBase.attr(
+ visible=false,
+ scaleanchor="y",
+ range=[0, width]
),
- yaxis = PlotlyBase.attr(
- visible = false,
- range = [-height,0]
+ yaxis=PlotlyBase.attr(
+ visible=false,
+ range=[-height, 0]
),
- margin = attr(l = 0, r = 0, t = 0, b = 0, pad = 0)
+ margin=attr(l=0, r=0, t=0, b=0, pad=0)
)
# Create the trace for the main image
trace = PlotlyBase.heatmap(
- z = elevation,
- x = X,
- y = -Y,
- name = "",
- hoverinfo = "x+y",
- showlegend = false,
- colorscale = "Viridis",
- showscale = false
+ z=elevation,
+ x=X,
+ y=-Y,
+ name="",
+ hoverinfo="x+y",
+ showlegend=false,
+ colorscale="Viridis",
+ showscale=false
)
plotdata = [trace]
@@ -188,10 +247,18 @@ function loadImgPlot(interfaceImg::String, overlayImg::String, imgTrans::Float64
return plotdata, plotlayout, width, height
end
-# loadContourPlot recieves the local directory of the image as a string,
-# returns the layout and data for the contour plotly plot
-# this function loads the image and applies a gaussian filter
-# to smoothen it and loads it into a plot
+"""
+ loadContourPlot(interfaceImg::String)
+
+Loads an image, smooths it, and creates a Plotly contour plot.
+
+# Arguments
+- `interfaceImg`: Path to the image file.
+
+# Returns
+- `plotdata`: A vector containing the Plotly contour trace.
+- `plotlayout`: The Plotly layout for the plot.
+"""
function loadContourPlot(interfaceImg::String)
# Load the image
cleaned_img=replace(interfaceImg, r"\?.*" => "")
@@ -258,10 +325,18 @@ function loadContourPlot(interfaceImg::String)
return plotdata, plotlayout
end
-# loadSurfacePlot recieves the local directory of the image as a string,
-# returns the layout and data for the surface plotly plot
-# this function loads the image and applies a gaussian filter
-# to smoothen it and loads it into a 3D plot
+"""
+ loadSurfacePlot(interfaceImg::String)
+
+Loads an image, smooths it, and creates a 3D Plotly surface plot.
+
+# Arguments
+- `interfaceImg`: Path to the image file.
+
+# Returns
+- `plotdata`: A vector containing the Plotly surface trace.
+- `plotlayout`: The Plotly layout for the 3D plot.
+"""
function loadSurfacePlot(interfaceImg::String)
# Load the image
cleaned_img=replace(interfaceImg, r"\?.*" => "")
@@ -347,28 +422,51 @@ function loadSurfacePlot(interfaceImg::String)
return plotdata, plotlayout
end
-# This function recieves the x and y coords currently selected, and the dimentions of
-# the image to create two traces that will display in a cross section
+"""
+ crossLinesPlot(x, y, maxwidth, maxheight)
+
+Creates two line traces for a crosshair indicator on a plot.
+
+# Arguments
+- `x`: The x-coordinate of the crosshair center.
+- `y`: The y-coordinate of the crosshair center.
+- `maxwidth`: The width of the plot area.
+- `maxheight`: The height of the plot area.
+
+# Returns
+- `trace1`: The horizontal line trace.
+- `trace2`: The vertical line trace.
+"""
function crossLinesPlot(x, y, maxwidth, maxheight)
# Define the coordinates for the two lines
- l1_x=[0, maxwidth]
- l1_y=[y, y]
- l2_x=[x, x]
- l2_y=[0, maxheight]
+ l1_x = [0, maxwidth]
+ l1_y = [y, y]
+ l2_x = [x, x]
+ l2_y = [0, maxheight]
# Create the line traces
- trace1=PlotlyBase.scatter(x=l1_x, y=l1_y, mode="lines",line=attr(color="red", width=0.5),name="Line X",showlegend=false)
- trace2=PlotlyBase.scatter(x=l2_x, y=l2_y, mode="lines",line=attr(color="red", width=0.5),name="Line Y",showlegend=false)
+ trace1 = PlotlyBase.scatter(x=l1_x, y=l1_y, mode="lines", line=attr(color="red", width=0.5), name="Line X", showlegend=false)
+ trace2 = PlotlyBase.scatter(x=l2_x, y=l2_y, mode="lines", line=attr(color="red", width=0.5), name="Line Y", showlegend=false)
return trace1, trace2
end
-# This function is used for giving colorbar values a visual format
-# that shortens long values giving them scientific notation
+"""
+ log_tick_formatter(values::Vector{Float64})
+
+Formats a vector of numbers into strings with a custom scientific notation for use as tick labels.
+For example, 1000 becomes "100x10¹" and 0.01 becomes "1.0x10⁻²".
+
+# Arguments
+- `values`: A vector of `Float64` values to format.
+
+# Returns
+- A vector of formatted strings.
+"""
function log_tick_formatter(values::Vector{Float64})
# Initialize exponents dictionary
- exponents=zeros(Int, length(values))
- formValues=zeros(Float64, length(values))
+ exponents = zeros(Int, length(values))
+ formValues = zeros(Float64, length(values))
for i in 1:length(values)
value = values[i]
if value >= 1000 # positive formatting for notation
@@ -382,78 +480,36 @@ function log_tick_formatter(values::Vector{Float64})
exponents[i] -= 1
end
end
- formValues[i]=value
+ formValues[i] = value
end
return map((v, e) -> e == 0 ? "$(round(v, sigdigits=2))" : "$(round(v, sigdigits=2))x10" * Makie.UnicodeFun.to_superscript(e), formValues, exponents)
end
-function generate_colorbar_image(slice_data::AbstractMatrix, color_levels::Int, output_path::String; use_triq::Bool=false, triq_prob::Float64=0.98)
- # 1. Determine bounds based on whether TrIQ is used
- min_val, max_val = if use_triq
- MSI_src.get_outlier_thres(slice_data, triq_prob)
+"""
+ meanSpectrumPlot(data::MSIData, dataset_name::String="")
+
+Generates a plot of the mean spectrum from MSIData.
+
+# Arguments
+- `data`: The `MSIData` object.
+- `dataset_name`: Optional name of the dataset for the plot title.
+
+# Returns
+- `plotdata`: A vector containing the Plotly trace.
+- `plotlayout`: The Plotly layout for the plot.
+- `xSpectraMz`: The m/z values of the spectrum.
+- `ySpectraMz`: The intensity values of the spectrum.
+"""
+function meanSpectrumPlot(data::MSIData, dataset_name::String=""; mask_path::Union{String, Nothing}=nothing)
+ # Determine base title based on mask usage
+ base_title = if mask_path !== nothing
+ "Masked Average Spectrum"
else
- extrema(slice_data)
- end
-
- # 2. Replicate the tick calculation logic from plot_slices
- bins = color_levels
- 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
+ "Average Spectrum"
end
- exponent = level_range > 0 ? floor(log10(level_range)) / 3 : 0
- scale = 10^(3 * exponent)
- scaled_levels = levels ./ scale
+ title_text = isempty(dataset_name) ? base_title : "$base_title for: $dataset_name"
- 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]
-
- # 3. Create and save the colorbar image
- fig = Figure(size=(150, 250))
- Colorbar(fig[1, 1],
- colormap=cgrad(:viridis, bins, categorical=true),
- # limits=(min_val, max_val),
- limits=(levels[1], levels[end]),
- label=(scale == 1 ? "Intensity" : "Intensity ×10^$(round(Int, 3 * exponent))"),
- ticks=(tick_positions, tick_labels),
- labelsize=20,
- ticklabelsize=16
- )
- save(output_path, fig)
-end
-
-# meanSpectrumPlot recieves the local directory of the image as a string,
-# returns the layout and data for the surface plotly plot
-# this function loads the spectra data and makes a mean to display
-# its values in the spectrum plot
-function meanSpectrumPlot(data::MSIData, dataset_name::String="")
- title_text = isempty(dataset_name) ? "Average Spectrum Plot" : "Average Spectrum for: $dataset_name"
layout = PlotlyBase.Layout(
title=PlotlyBase.attr(
text=title_text,
@@ -477,42 +533,76 @@ function meanSpectrumPlot(data::MSIData, dataset_name::String="")
)
# Use the new, efficient function from the backend
- xSpectraMz, ySpectraMz = get_average_spectrum(data)
+ xSpectraMz, ySpectraMz = get_average_spectrum(data, mask_path=mask_path)
- if isempty(xSpectraMz)
+ if isempty(xSpectraMz) || isempty(ySpectraMz)
@warn "Average spectrum is empty."
- trace = PlotlyBase.scatter(x=Float64[], y=Float64[])
+ trace = PlotlyBase.stem(x=Float64[], y=Float64[])
+ # Update title to indicate empty spectrum
+ layout.title.text = "Empty " * layout.title.text
else
- trace = PlotlyBase.scatter(x=xSpectraMz, y=ySpectraMz, marker=attr(size=1, color="blue", opacity=0.5), name="Average", hoverinfo="x",hovertemplate="
m/z: %{x:.4f}
")
+ trace = PlotlyBase.stem(x=xSpectraMz, y=ySpectraMz, marker=attr(size=1, color="blue", opacity=0.5), name="Average", hoverinfo="x", hovertemplate="
m/z: %{x:.4f}
")
end
-
plotdata = [trace]
plotlayout = layout
return plotdata, plotlayout, xSpectraMz, ySpectraMz
end
-function xySpectrumPlot(data::MSIData, xCoord::Int, yCoord::Int, imgWidth::Int, imgHeight::Int, dataset_name::String="")
+"""
+ xySpectrumPlot(data::MSIData, xCoord::Int, yCoord::Int, imgWidth::Int, imgHeight::Int, dataset_name::String="")
+
+Generates a plot for the spectrum at a specific coordinate (for imaging data) or index (for non-imaging data).
+
+# Arguments
+- `data`: The `MSIData` object.
+- `xCoord`: The x-coordinate or spectrum index.
+- `yCoord`: The y-coordinate (used for imaging data).
+- `imgWidth`: The width of the MSI image.
+- `imgHeight`: The height of the MSI image.
+- `dataset_name`: Optional name of the dataset for the plot title.
+
+# Returns
+- `plotdata`: A vector containing the Plotly trace.
+- `plotlayout`: The Plotly layout for the plot.
+- `mz`: The m/z values of the spectrum.
+- `intensity`: The intensity values of the spectrum.
+"""
+function xySpectrumPlot(data::MSIData, xCoord::Int, yCoord::Int, imgWidth::Int, imgHeight::Int, dataset_name::String=""; mask_path::Union{String, Nothing}=nothing)
local mz::AbstractVector, intensity::AbstractVector
local plot_title::String
is_imaging = data.source isa ImzMLSource
if is_imaging
- # For imaging data, use (X, Y) coordinates
x = clamp(xCoord, 1, imgWidth)
y = clamp(yCoord, 1, imgHeight)
-
- process_spectrum(data, Int(x), Int(y)) do recieved_mz, recieved_intensity
- mz = recieved_mz
- intensity = recieved_intensity
+
+ if mask_path !== nothing
+ mask_matrix = load_and_prepare_mask(mask_path, (imgWidth, imgHeight))
+ if !mask_matrix[y, x]
+ @warn "Coordinate ($x, $y) is outside the specified mask."
+ mz, intensity = Float64[], Float64[]
+ base_title = "Spectrum Outside Mask at ($x, $y)"
+ else
+ process_spectrum(data, Int(x), Int(y)) do recieved_mz, recieved_intensity
+ mz = recieved_mz
+ intensity = recieved_intensity
+ end
+ base_title = "Masked Spectrum at ($x, $y)"
+ end
+ else
+ process_spectrum(data, Int(x), Int(y)) do recieved_mz, recieved_intensity
+ mz = recieved_mz
+ intensity = recieved_intensity
+ end
+ base_title = "Spectrum at ($x, $y)"
end
- base_title = "Spectrum at ($x, $y)"
else
# For non-imaging data, treat xCoord as the spectrum index
index = clamp(xCoord, 1, length(data.spectra_metadata))
-
- process_spectrum(data, index) do recieved_mz, recieved_intensity
+
+ process_spectrum(data, index) do recieved_mz, recieved_intensity
mz = recieved_mz
intensity = recieved_intensity
end
@@ -542,20 +632,43 @@ function xySpectrumPlot(data::MSIData, xCoord::Int, yCoord::Int, imgWidth::Int,
),
margin=attr(l=0, r=0, t=120, b=0, pad=0)
)
-
+
# Downsample for plotting performance
mz_down, int_down = MSI_src.downsample_spectrum(mz, intensity)
- trace = PlotlyBase.scatter(x=mz_down, y=int_down, marker=attr(size=1, color="blue", opacity=0.5), name="Spectrum", hoverinfo="x", hovertemplate="
m/z: %{x:.4f}
")
-
+ trace = PlotlyBase.stem(x=mz_down, y=int_down, marker=attr(size=1, color="blue", opacity=0.5), name="Spectrum", hoverinfo="x", hovertemplate="
m/z: %{x:.4f}
")
+
plotdata = [trace]
plotlayout = layout
-
+
return plotdata, plotlayout, mz, intensity
end
-function sumSpectrumPlot(data::MSIData, dataset_name::String="")
- title_text = isempty(dataset_name) ? "Total Spectrum Plot" : "Total Spectrum for: $dataset_name"
+"""
+ sumSpectrumPlot(data::MSIData, dataset_name::String="")
+
+Generates a plot of the total (summed) spectrum from MSIData.
+
+# Arguments
+- `data`: The `MSIData` object.
+- `dataset_name`: Optional name of the dataset for the plot title.
+
+# Returns
+- `plotdata`: A vector containing the Plotly trace.
+- `plotlayout`: The Plotly layout for the plot.
+- `xSpectraMz`: The m/z values of the spectrum.
+- `ySpectraMz`: The intensity values of the spectrum.
+"""
+function sumSpectrumPlot(data::MSIData, dataset_name::String=""; mask_path::Union{String, Nothing}=nothing)
+ # Determine base title based on mask usage
+ base_title = if mask_path !== nothing
+ "Masked Total Spectrum"
+ else
+ "Total Spectrum"
+ end
+
+ title_text = isempty(dataset_name) ? base_title : "$base_title for: $dataset_name"
+
layout = PlotlyBase.Layout(
title=PlotlyBase.attr(
text=title_text,
@@ -579,13 +692,15 @@ function sumSpectrumPlot(data::MSIData, dataset_name::String="")
)
# Use the get_total_spectrum function from the backend
- xSpectraMz, ySpectraMz = get_total_spectrum(data)
+ xSpectraMz, ySpectraMz, num_spectra = get_total_spectrum(data, mask_path=mask_path)
- if isempty(xSpectraMz)
+ if isempty(xSpectraMz) || isempty(ySpectraMz)
@warn "Total spectrum is empty."
- trace = PlotlyBase.scatter(x=Float64[], y=Float64[])
+ trace = PlotlyBase.stem(x=Float64[], y=Float64[])
+ # Update title to indicate empty spectrum
+ layout.title.text = "Empty " * layout.title.text
else
- trace = PlotlyBase.scatter(x=xSpectraMz, y=ySpectraMz, marker=attr(size=1, color="blue", opacity=0.5), name="Total", hoverinfo="x",hovertemplate="
m/z: %{x:.4f}
")
+ trace = PlotlyBase.stem(x=xSpectraMz, y=ySpectraMz, marker=attr(size=1, color="blue", opacity=0.5), name="Total", hoverinfo="x", hovertemplate="
m/z: %{x:.4f}
")
end
plotdata = [trace]
@@ -593,30 +708,32 @@ function sumSpectrumPlot(data::MSIData, dataset_name::String="")
return plotdata, plotlayout, xSpectraMz, ySpectraMz
end
+"""
+ warmup_init()
+
+Performs pre-compilation of key functions at application startup to reduce first-use latency.
+This is run asynchronously and should not block application startup.
+"""
function warmup_init()
@async begin
println("Pre-compiling functions at startup...")
-
- # Create a dummy MSIData object to be used for pre-compilation
- # dummy_source = ImzMLSource("dummy.ibd", Float32, Float32)
- # dummy_meta = MSI_src.SpectrumMetadata(0,0,"",MSI_src.UNKNOWN, MSI_src.SpectrumAsset(Float32,false,0,0,:mz), MSI_src.SpectrumAsset(Float32,false,0,0,:intensity))
- # dummy_msi_data = MSIData(dummy_source, [dummy_meta], (1,1), zeros(Int,1,1), 0)
- # Pre-compile functions from btnSearch
- # try OpenMSIData("dummy.imzML") catch end
- # try precompute_analytics(dummy_msi_data) catch end
+ # Pre-compile image processing and plotting functions
+ try
+ TrIQ(zeros(10, 10), 256, 0.98)
+ catch
+ end
+ try
+ quantize_intensity(zeros(10, 10), 256)
+ catch
+ end
- # Pre-compile functions from mainProcess
- # try get_mz_slice(dummy_msi_data, 1.0, 1.0) catch end
- try TrIQ(zeros(10,10), 256, 0.98) catch end
- try quantize_intensity(zeros(10,10), 256) catch end
-
dummy_bmp_path = joinpath("public", "dummy.bmp")
dummy_png_path = joinpath("public", "dummy.png")
- try
+ try
save_bitmap(dummy_bmp_path, zeros(UInt8, 10, 10), ViridisPalette)
loadImgPlot("/dummy.bmp")
- generate_colorbar_image(zeros(10,10), 256, dummy_png_path)
+ generate_colorbar_image(zeros(10, 10), 256, dummy_png_path)
catch e
@warn "Pre-compilation step failed (this is expected if dummy files can't be created/read)"
finally
@@ -628,18 +745,41 @@ function warmup_init()
end
end
+"""
+ load_registry(registry_path)
+
+Loads the dataset registry from a JSON file.
+
+# Arguments
+- `registry_path`: Path to the `registry.json` file.
+
+# Returns
+- A dictionary containing the registry data. Returns an empty dictionary if the file doesn't exist or fails to parse.
+"""
function load_registry(registry_path)
if isfile(registry_path)
try
- return JSON.parsefile(registry_path, dicttype=Dict{String, Any})
+ return JSON.parsefile(registry_path, dicttype=Dict{String,Any})
catch e
@error "Failed to parse registry.json: $e"
- return Dict{String, Any}()
+ return Dict{String,Any}()
end
end
- return Dict{String, Any}()
+ return Dict{String,Any}()
end
+"""
+ extract_metadata(msi_data::MSIData, source_path::String)
+
+Extracts key metadata from an `MSIData` object for display.
+
+# Arguments
+- `msi_data`: The `MSIData` object.
+- `source_path`: The path to the source data file.
+
+# Returns
+- A dictionary containing summary statistics and metadata.
+"""
function extract_metadata(msi_data::MSIData, source_path::String)
df = msi_data.spectrum_stats_df
if df === nothing
@@ -671,7 +811,7 @@ function extract_metadata(msi_data::MSIData, source_path::String)
centroid_count = count(==(MSI_src.CENTROID), df.Mode)
profile_count = count(==(MSI_src.PROFILE), df.Mode)
unknown_count = count(==(MSI_src.UNKNOWN), df.Mode)
-
+
push!(summary_stats, Dict("parameter" => "Centroid Spectra", "value" => string(centroid_count)))
push!(summary_stats, Dict("parameter" => "Profile Spectra", "value" => string(profile_count)))
if unknown_count > 0
@@ -686,18 +826,39 @@ function extract_metadata(msi_data::MSIData, source_path::String)
)
end
+"""
+ update_registry(registry_path, dataset_name, source_path, metadata=nothing, is_imzML=false)
+
+Adds or updates an entry in the dataset registry JSON file.
+
+# Arguments
+- `registry_path`: Path to the `registry.json` file.
+- `dataset_name`: The name of the dataset.
+- `source_path`: The path to the source data file.
+- `metadata`: Optional dictionary of metadata to store.
+- `is_imzML`: Boolean indicating if the source is an imzML file.
+"""
function update_registry(registry_path, dataset_name, source_path, metadata=nothing, is_imzML=false)
registry = load_registry(registry_path)
- entry = Dict{String, Any}( # Explicitly type the dictionary to allow mixed value types
- "source_path" => source_path,
- "processed_date" => string(now()),
- "is_imzML" => is_imzML
- )
+
+ # Get existing entry if it exists, otherwise create new one
+ existing_entry = get(registry, dataset_name, Dict{String,Any}())
+
+ # Start with existing data and update only the basic fields
+ entry = copy(existing_entry)
+ entry["source_path"] = source_path
+ entry["processed_date"] = string(now())
+ entry["is_imzML"] = is_imzML
+
+ # Only update metadata if provided
if metadata !== nothing
entry["metadata"] = metadata
end
- registry[dataset_name] = entry
+ # Note: has_mask and mask_path are preserved from existing_entry if they exist
+
+ registry[dataset_name] = entry
+
try
open(registry_path, "w") do f
JSON.print(f, registry, 4)
@@ -707,12 +868,30 @@ function update_registry(registry_path, dataset_name, source_path, metadata=noth
end
end
-function process_file_safely(file_path, masses, params, progress_message_ref, overall_progress_ref)
+"""
+ process_file_safely(file_path, masses, params, progress_message_ref, overall_progress_ref)
+
+Safely processes a single MSI data file, generating and saving m/z slices.
+
+This function handles loading data, generating slices, saving results as bitmaps,
+and updating the data registry. It includes error handling and memory cleanup.
+
+# Arguments
+- `file_path`: Path to the `.imzML` file.
+- `masses`: A vector of m/z values to generate slices for.
+- `params`: A structure or dictionary containing processing parameters.
+- `progress_message_ref`: A reference to update with progress messages.
+- `overall_progress_ref`: A reference to update with overall progress.
+
+# Returns
+- A tuple `(success::Bool, message::String)`.
+"""
+function process_file_safely(file_path, masses, params, progress_message_ref, overall_progress_ref; use_mask::Bool=false)
local_msi_data = nothing
dataset_name = replace(basename(file_path), r"\.imzML$"i => "")
output_dir = joinpath("public", dataset_name)
println("Processing: $dataset_name -> $output_dir")
-
+
try
# --- Load Data ---
progress_message_ref = "Loading: $(basename(file_path))"
@@ -722,19 +901,43 @@ function process_file_safely(file_path, masses, params, progress_message_ref, ov
return (false, "Skipped: Not an imzML file")
end
- # --- Generate Slices (this will call precompute_analytics if needed) ---
- progress_message_ref = "Generating $(length(masses)) slices for $(dataset_name)..."
- slice_dict = get_multiple_mz_slices(local_msi_data, masses, params.tolerance)
+ # --- Get mask path and load mask_matrix only if use_mask is true ---
+ local mask_matrix_for_triq::Union{BitMatrix, Nothing} = nothing
+ mask_path = nothing
+ if use_mask
+ registry = load_registry(params.registry)
+ entry = get(registry, dataset_name, nothing)
+ if entry !== nothing && get(entry, "has_mask", false)
+ mask_path_candidate = get(entry, "mask_path", "")
+ if isfile(mask_path_candidate)
+ mask_path = mask_path_candidate
+ # Load the mask matrix here
+ mask_matrix_for_triq = load_and_prepare_mask(mask_path, (local_msi_data.image_dims[1], local_msi_data.image_dims[2]))
+ println("DEBUG: Using mask: $(mask_path) with dimensions $(size(mask_matrix_for_triq))")
+ else
+ @warn "Mask enabled but file not found: $(mask_path_candidate). Clearing invalid mask entry."
+ # Clear invalid mask entry
+ update_registry_mask_fields(params.registry, dataset_name, false, "")
+ mask_path = nothing
+ end
+ else
+ @warn "Mask enabled but no valid mask entry found for: $(dataset_name)"
+ end
+ end
- # --- Extract metadata *after* it has been computed ---
+ # --- Generate Slices ---
+ progress_message_ref = "Generating $(length(masses)) slices for $(dataset_name)..."
+ slice_dict = get_multiple_mz_slices(local_msi_data, masses, params.tolerance, mask_path=mask_path)
+
+ # --- Extract metadata ---
metadata = extract_metadata(local_msi_data, file_path)
# --- Save Slices ---
- mkpath(output_dir) # Ensure output directory exists
-
+ mkpath(output_dir)
+
for (mass_idx, mass) in enumerate(masses)
progress_message_ref = "File $(params.fileIdx)/$(params.nFiles): Saving slice for m/z=$mass"
-
+
slice = slice_dict[mass]
text_nmass = replace(string(mass), "." => "_")
bitmap_filename = params.triqE ? "TrIQ_$(text_nmass).bmp" : "MSI_$(text_nmass).bmp"
@@ -744,7 +947,7 @@ function process_file_safely(file_path, masses, params, progress_message_ref, ov
sliceQuant = zeros(UInt8, size(slice))
@warn "No intensity data for m/z = $mass in $(dataset_name)"
else
- sliceQuant = params.triqE ? TrIQ(slice, params.colorL, params.triqP) : quantize_intensity(slice, params.colorL)
+ sliceQuant = params.triqE ? TrIQ(slice, params.colorL, params.triqP, mask_matrix=mask_matrix_for_triq) : quantize_intensity(slice, params.colorL, mask_matrix=mask_matrix_for_triq)
if params.medianF
sliceQuant = round.(UInt8, median_filter(sliceQuant))
end
@@ -752,10 +955,11 @@ function process_file_safely(file_path, masses, params, progress_message_ref, ov
save_bitmap(joinpath(output_dir, bitmap_filename), sliceQuant, ViridisPalette)
if !all(iszero, slice)
- generate_colorbar_image(slice, params.colorL, joinpath(output_dir, colorbar_filename); use_triq=params.triqE, triq_prob=params.triqP)
+ generate_colorbar_image(slice, params.colorL, joinpath(output_dir, colorbar_filename);
+ use_triq=params.triqE, triq_prob=params.triqP, mask_path=mask_path)
end
end
-
+
is_imzML = local_msi_data.source isa ImzMLSource
update_registry(params.registry, dataset_name, file_path, metadata, is_imzML)
return (true, "")
@@ -773,4 +977,150 @@ function process_file_safely(file_path, masses, params, progress_message_ref, ov
ccall(:malloc_trim, Int32, (Int32,), 0)
end
end
+end
+
+function update_registry_mask_fields(registry_path, dataset_name, has_mask, mask_path)
+ registry = load_registry(registry_path)
+
+ if haskey(registry, dataset_name)
+ registry[dataset_name]["has_mask"] = has_mask
+ registry[dataset_name]["mask_path"] = mask_path
+ else
+ # Create new entry if it doesn't exist
+ registry[dataset_name] = Dict{String,Any}(
+ "has_mask" => has_mask,
+ "mask_path" => mask_path,
+ "source_path" => "",
+ "processed_date" => string(now()),
+ "is_imzML" => false
+ )
+ end
+
+ try
+ open(registry_path, "w") do f
+ JSON.print(f, registry, 4)
+ end
+ catch e
+ @error "Failed to write to registry.json: $e"
+ end
+end
+
+"""
+ loadSurfacePlot(interfaceImg::String, mask_path::String)
+
+Loads an image, smooths it, applies a mask, crops the data to the masked region, and creates a 3D Plotly surface plot.
+
+# Arguments
+- `interfaceImg`: Path to the image file.
+- `mask_path`: Path to a mask file. Areas outside the mask will be removed and the plot axes will be cropped to the masked region.
+
+# Returns
+- `plotdata`: A vector containing the Plotly surface trace.
+- `plotlayout`: The Plotly layout for the 3D plot.
+"""
+function loadSurfacePlot(interfaceImg::String, mask_path::String)
+ # Load the image
+ cleaned_img = replace(interfaceImg, r"\?.*" => "")
+ cleaned_img = lstrip(cleaned_img, '/')
+ var = joinpath("./public", cleaned_img)
+ img = load(var)
+ img_gray = Gray.(img) # Convert to grayscale
+ img_array = Array(img_gray)
+ elevation = Float32.(Array(img_array)) ./ 255.0 # Normalize between 0 and 1
+
+ # Smooth the image
+ sigma = 3.0
+ kernel = Kernel.gaussian(sigma)
+ elevation_smoothed = imfilter(elevation, kernel)
+
+ # --- APPLY MASK and CROP DATA ---
+ mask_matrix = load_and_prepare_mask(mask_path, (size(elevation_smoothed, 2), size(elevation_smoothed, 1)))
+ elevation_smoothed[.!mask_matrix] .= NaN32
+
+ non_nan_indices = findall(!isnan, elevation_smoothed)
+ if isempty(non_nan_indices)
+ # Return an empty plot if mask removes everything
+ return [PlotlyBase.surface()], PlotlyBase.Layout(title="Empty plot: Mask covered all data")
+ end
+
+ row_indices = [idx[1] for idx in non_nan_indices]
+ col_indices = [idx[2] for idx in non_nan_indices]
+
+ y_range = minimum(row_indices):maximum(row_indices)
+ x_range = minimum(col_indices):maximum(col_indices)
+
+ cropped_elevation = elevation_smoothed[y_range, x_range]
+ # ---
+
+ # --- DOWNSAMPLE CROPPED DATA ---
+ cropped_elevation = downsample_image(cropped_elevation, 256)
+ # ---
+
+ # Create the X, Y meshgrid coordinates for the CROPPED data
+ x_coords = x_range
+ y_coords = y_range
+ X = repeat(reshape(x_coords, 1, length(x_coords)), length(y_coords), 1)
+ Y = repeat(reshape(y_coords, length(y_coords), 1), 1, length(x_coords))
+
+ # Define tick values and text for colorbars, ignoring NaNs
+ non_nan_values = filter(!isnan, cropped_elevation)
+ min_val = isempty(non_nan_values) ? 0.0 : minimum(non_nan_values)
+ max_val = isempty(non_nan_values) ? 1.0 : maximum(non_nan_values)
+ tickV = range(min_val, stop=max_val, length=8)
+ tickT = log_tick_formatter(collect(tickV))
+
+ # Transpose the elevation_smoothed array if Y axis is longer than X axis to fix chopping
+ cropped_elevation = transpose(cropped_elevation)
+ if size(cropped_elevation, 1) < size(cropped_elevation, 2)
+ Y = -Y
+ else
+ X = -X
+ end
+
+ # Calculate the number of ticks and aspect ratio for the 3d plot
+ x_nticks = min(20, length(x_coords))
+ y_nticks = min(20, length(y_coords))
+ z_nticks = 5
+ aspect_ratio = attr(x=1, y=length(y_coords) / length(x_coords), z=0.5)
+
+ # Define the layout for the 3D plot
+ layout3D = PlotlyBase.Layout(
+ title=PlotlyBase.attr(
+ text="3D surface plot of $cleaned_img (masked, cropped, downsampled)",
+ font=PlotlyBase.attr(
+ family="Roboto, Lato, sans-serif",
+ size=18,
+ color="black"
+ )
+ ),
+ scene=attr(
+ xaxis_nticks=x_nticks,
+ yaxis_nticks=y_nticks,
+ zaxis_nticks=z_nticks,
+ camera=attr(eye=attr(x=0, y=1, z=0.5)),
+ aspectratio=aspect_ratio
+ ),
+ margin=attr(l=0, r=0, t=120, b=0, pad=0)
+ )
+
+ trace3D = PlotlyBase.surface(
+ x=X[1, :],
+ y=Y[:, 1],
+ z=cropped_elevation,
+ contours_z=attr(
+ show=true,
+ usecolormap=true,
+ highlightcolor="limegreen",
+ project_z=true
+ ),
+ colorscale="Viridis",
+ colorbar=attr(
+ tickvals=tickV,
+ ticktext=tickT,
+ nticks=8
+ )
+ )
+ plotdata = [trace3D]
+ plotlayout = layout3D
+ return plotdata, plotlayout
end
\ No newline at end of file
diff --git a/mask.jl b/mask.jl
index 349a5f2..1682805 100644
--- a/mask.jl
+++ b/mask.jl
@@ -9,14 +9,12 @@ using Statistics, NaturalSort, LinearAlgebra, StipplePlotly
using Base.Filesystem: mv
using MSI_src
-using .MSI_src: MSIData
+using .MSI_src: MSIData, process_image_pipeline
# Plot Handling
include("./julia_imzML_visual.jl")
# Image Processing Pipeline
-include("src/ImageProcessing.jl")
-using .ImageProcessing
using ImageBinarization
function load_and_binarize_mask(path)
@@ -46,7 +44,7 @@ function alter_image(img_path, otsu_scale, noise_size_percent, hole_size_percent
gray_img = Float32.(ensure_grayscale(original_img))
binary, noise_removed, holes_filled, smoothed =
- ImageProcessing.process_image_pipeline(gray_img;
+ process_image_pipeline(gray_img;
otsu_scale=otsu_scale, noise_size_percent=noise_size_percent,
hole_size_percent=hole_size_percent, smoothing=smoothing_level)
@@ -250,7 +248,6 @@ end
is_browsing_slices = true
is_editing_mask = false
folder_path = joinpath("public", selected_folder_main)
- println("Selected folder: $selected_folder_main")
if !isdir(folder_path)
imgInt = ""
msgimg = "Folder not found."
@@ -674,15 +671,17 @@ end
# Update registry directly in the handler
reg_path = abspath(joinpath(@__DIR__, "public", "registry.json"))
registry = isfile(reg_path) ? JSON.parsefile(reg_path) : Dict{String, Any}()
+
+ full_mask_path = abspath(final_path)
# Update the registry entry
if haskey(registry, selected_folder_main)
- registry[selected_folder_main]["mask_path"] = "/css/masks/$(final_mask_name)"
+ registry[selected_folder_main]["mask_path"] = full_mask_path
registry[selected_folder_main]["has_mask"] = true
else
# Create a new entry if folder doesn't exist in registry
registry[selected_folder_main] = Dict(
- "mask_path" => "/css/masks/$(final_mask_name)",
+ "mask_path" => full_mask_path,
"has_mask" => true,
"is_imzML" => true,
"processed_date" => string(Dates.now())
diff --git a/mask.jl.html b/mask.jl.html
index 5d5492d..60d7011 100644
--- a/mask.jl.html
+++ b/mask.jl.html
@@ -1,72 +1,91 @@
-
-
-
-
Mask Generation Workflow
-
-
-
Step 1: Select Slice
-
-
-
-
+
+
+
+
Mask Generation Workflow
-
-
Step 2: Create/Upload Mask
-
-
-
-
-
+
+
Step 1: Select Slice
+
+
+
+
-
-
Step 3: Adjust & Verify
-
-
Adjust automated processing for the initial mask:
-
-
-
-
-
-
-
Adjust slice overlay transparency:
-
-
+
+
Step 2: Create/Upload Mask
+
+
+
+
+
-
{{ mask_editor_message }}
-
-
-
+
+
Step 3: Adjust & Verify
+
+
Adjust automated processing for the initial mask:
+
+
+
+
+
+
+
Adjust slice overlay transparency:
+
+
+
+
{{ mask_editor_message }}
+
+
+
+
-
Mask Preview
-
-
+
+
-
+
Please select a dataset and provide a mask input to begin.
Use the controls on the left to load a slice or upload a mask image.
-
+
@@ -87,54 +106,59 @@
-
-
-
-
+
+
+
+
-
+
-
-
-
+
+
+
-
-
-
+
+
+
-
-
+
+
-
+
-
-
+ }" v-on:mousedown.prevent="startDrawing" v-on:mousemove.prevent="draw"
+ v-on:mouseup.prevent="stopDrawing()" v-on:mouseleave.prevent="stopDrawing()">
+
+
@@ -142,4 +166,4 @@
-->
-
+
\ No newline at end of file
diff --git a/src/ImageProcessing.jl b/src/ImageProcessing.jl
index 2f0f8a9..d47ea6a 100644
--- a/src/ImageProcessing.jl
+++ b/src/ImageProcessing.jl
@@ -1,16 +1,80 @@
-module ImageProcessing
+# src/ImageProcessing
using Images
using ImageBinarization
using ImageMorphology
using ImageComponentAnalysis
+using Colors # For converting to grayscale
export process_image_pipeline
+export load_and_prepare_mask # Export the new function
+
+"""
+ load_and_prepare_mask(mask_path::String, target_dims::Tuple{Int, Int})
+
+Loads a PNG image mask, converts it to a binary (Boolean) matrix,
+and resizes it to the specified `target_dims`. White pixels in the mask
+are considered `true` (part of the ROI), and black pixels are `false`.
+
+# Arguments
+- `mask_path`: Absolute path to the PNG mask file.
+- `target_dims`: A tuple `(width, height)` representing the desired output dimensions.
+
+# Returns
+- A `BitMatrix` of `target_dims` where `true` indicates the ROI.
+"""
+function load_and_prepare_mask(mask_path::String, target_dims::Tuple{Int, Int})
+ if !isfile(mask_path)
+ error("Mask file not found: $(mask_path)")
+ end
+
+ # Load the image
+ img = Images.load(mask_path)
+
+ # Convert to grayscale if it's a color image
+ gray_img = Gray.(img)
+
+ # Binarize using Otsu's method - white regions become 'true'
+ binary_img = binarize(gray_img, Otsu())
+
+ # Resize to target dimensions
+ resized_img = imresize(binary_img, (target_dims[2], target_dims[1]))
+
+ # Ensure the output is a BitMatrix, as imresize can change the type
+ return resized_img .> 0.5
+end
# ===================================================================
# CORE PROCESSING PIPELINE
# ===================================================================
+"""
+ process_image_pipeline(gray_img; otsu_scale=1.0, noise_size_percent=0.1, hole_size_percent=0.05, smoothing=2)
+
+Applies a multi-step image processing pipeline to a grayscale image to segment regions of interest.
+
+The pipeline consists of:
+1. **Binarization**: An adjusted Otsu's threshold is used to create a binary image.
+2. **Noise Removal**: Small white regions (noise) are removed using an area opening operation.
+3. **Hole Filling**: Small black regions (holes) within larger objects are filled.
+4. **Edge Smoothing**: The edges of the final regions are smoothed using a morphological closing operation.
+
+# Arguments
+- `gray_img`: The input grayscale image (`Matrix{<:Gray}`).
+
+# Keyword Arguments
+- `otsu_scale`: A factor to scale the automatically determined Otsu threshold. Values > 1.0 make the threshold stricter (less white), < 1.0 make it more lenient (more white). Default: `1.0`.
+- `noise_size_percent`: The percentage of the total image area used as a threshold to remove small white noise components. Default: `0.1`.
+- `hole_size_percent`: The percentage of the total image area used as a threshold to fill black holes in white components. Default: `0.05`.
+- `smoothing`: The size of the kernel for the final edge smoothing (closing) operation. Default: `2`.
+
+# Returns
+- A tuple containing four images representing the intermediate steps of the pipeline:
+ 1. `binary_img`: The result of the initial binarization.
+ 2. `noise_removed_img`: The image after noise removal.
+ 3. `holes_filled_img`: The image after filling holes.
+ 4. `smoothed_img`: The final smoothed image.
+"""
function process_image_pipeline(gray_img;
otsu_scale=1.0,
noise_size_percent=0.1,
@@ -43,5 +107,3 @@ function process_image_pipeline(gray_img;
# --- Return all intermediate steps for visualization ---
return binary_img, noise_removed_img, holes_filled_img, smoothed_img
end
-
-end
diff --git a/src/MSIData.jl b/src/MSIData.jl
index 709b00f..b5bc5a7 100644
--- a/src/MSIData.jl
+++ b/src/MSIData.jl
@@ -143,6 +143,45 @@ mutable struct MSIData
end
end
+"""
+ get_masked_spectrum_indices(data::MSIData, mask_matrix::BitMatrix) -> Set{Int}
+
+Converts a 2D boolean mask matrix (e.g., loaded from a PNG) into a `Set` of linear
+spectrum indices that fall within the `true` regions of the mask.
+
+# Arguments
+- `data`: The `MSIData` object containing the `coordinate_map`.
+- `mask_matrix`: A `BitMatrix` where `true` indicates a pixel is part of the ROI.
+
+# Returns
+- A `Set{Int}` containing the linear indices of spectra within the mask.
+"""
+function get_masked_spectrum_indices(data::MSIData, mask_matrix::BitMatrix)
+ if data.coordinate_map === nothing
+ error("Coordinate map not available. Cannot apply mask to non-imaging data.")
+ end
+
+ width, height = data.image_dims
+ mask_height, mask_width = size(mask_matrix)
+
+ if mask_width != width || mask_height != height
+ error("Mask dimensions ($(mask_width)x$(mask_height)) do not match image dimensions ($(width)x$(height)).")
+ end
+
+ masked_indices = Set{Int}()
+ for y in 1:height
+ for x in 1:width
+ if mask_matrix[y, x] # If this pixel is part of the mask
+ idx = data.coordinate_map[x, y]
+ if idx != 0 # Ensure there's an actual spectrum at this coordinate
+ push!(masked_indices, idx)
+ end
+ end
+ end
+ end
+ return masked_indices
+end
+
# --- Internal function for reading binary data --- #
@@ -513,7 +552,7 @@ It uses a fast, two-pass approach:
This function is highly optimized for `.imzML` by leveraging direct binary
reading and optimized binning logic. It is called by `get_total_spectrum`.
"""
-function get_total_spectrum_imzml(msi_data::MSIData; num_bins::Int=2000)
+function get_total_spectrum_imzml(msi_data::MSIData; num_bins::Int=2000, masked_indices::Union{Set{Int}, Nothing}=nothing)
println("Calculating total spectrum for imzML (2-pass method)...")
total_start_time = time_ns()
@@ -529,7 +568,7 @@ function get_total_spectrum_imzml(msi_data::MSIData; num_bins::Int=2000)
println(" Pass 1: Finding global m/z range...")
g_min_mz = Inf
g_max_mz = -Inf
- _iterate_spectra_fast(msi_data) do idx, mz, _
+ _iterate_spectra_fast(msi_data, masked_indices === nothing ? nothing : collect(masked_indices)) do idx, mz, _
if !isempty(mz)
local_min, local_max = extrema(mz)
g_min_mz = min(g_min_mz, local_min)
@@ -539,7 +578,7 @@ function get_total_spectrum_imzml(msi_data::MSIData; num_bins::Int=2000)
pass1_duration = (time_ns() - pass1_start_time) / 1e9
if !isfinite(g_min_mz)
@warn "Could not determine a valid m/z range for imzML. All spectra might be empty."
- return (Float64[], Float64[])
+ return (Float64[], Float64[], 0)
end
# Use and cache the result
@@ -560,7 +599,9 @@ function get_total_spectrum_imzml(msi_data::MSIData; num_bins::Int=2000)
# 3. Second Pass: Optimized binning
pass2_start_time = time_ns()
println(" Pass 2: Summing intensities into $num_bins bins...")
- _iterate_spectra_fast(msi_data) do idx, mz, intensity
+ num_spectra_processed = 0
+ _iterate_spectra_fast(msi_data, masked_indices === nothing ? nothing : collect(masked_indices)) do idx, mz, intensity
+ num_spectra_processed += 1
if isempty(mz)
return
end
@@ -594,7 +635,7 @@ function get_total_spectrum_imzml(msi_data::MSIData; num_bins::Int=2000)
println("----------------------------------------\n")
println("Total spectrum calculation complete for imzML.")
- return (collect(mz_bins), intensity_sum)
+ return (collect(mz_bins), intensity_sum, num_spectra_processed)
end
"""
@@ -608,7 +649,7 @@ It uses a two-pass approach analogous to the `imzML` implementation:
This function is called by `get_total_spectrum`.
"""
-function get_total_spectrum_mzml(msi_data::MSIData; num_bins::Int=2000)
+function get_total_spectrum_mzml(msi_data::MSIData; num_bins::Int=2000, masked_indices::Union{Set{Int}, Nothing}=nothing)
println("Calculating total spectrum for mzML (2-pass method)...")
total_start_time = time_ns()
@@ -624,7 +665,7 @@ function get_total_spectrum_mzml(msi_data::MSIData; num_bins::Int=2000)
println(" Pass 1: Finding global m/z range...")
g_min_mz = Inf
g_max_mz = -Inf
- _iterate_spectra_fast(msi_data) do idx, mz, intensity
+ _iterate_spectra_fast(msi_data, masked_indices === nothing ? nothing : collect(masked_indices)) do idx, mz, intensity
if !isempty(mz)
local_min, local_max = extrema(mz)
g_min_mz = min(g_min_mz, local_min)
@@ -635,7 +676,7 @@ function get_total_spectrum_mzml(msi_data::MSIData; num_bins::Int=2000)
if !isfinite(g_min_mz)
@warn "Could not determine a valid m/z range for mzML. All spectra might be empty."
- return (Float64[], Float64[])
+ return (Float64[], Float64[], 0)
end
# Use and cache the result
@@ -658,7 +699,9 @@ function get_total_spectrum_mzml(msi_data::MSIData; num_bins::Int=2000)
inv_bin_step = 1.0 / bin_step # Precompute reciprocal
min_mz = global_min_mz
- _iterate_spectra_fast(msi_data) do idx, mz, intensity
+ num_spectra_processed = 0
+ _iterate_spectra_fast(msi_data, masked_indices === nothing ? nothing : collect(masked_indices)) do idx, mz, intensity
+ num_spectra_processed += 1
if isempty(mz)
return
end
@@ -687,7 +730,7 @@ function get_total_spectrum_mzml(msi_data::MSIData; num_bins::Int=2000)
println("----------------------------------------\n")
println("Total spectrum calculation complete for mzML.")
- return (collect(mz_bins), intensity_sum)
+ return (collect(mz_bins), intensity_sum, num_spectra_processed)
end
"""
@@ -699,11 +742,18 @@ This function dispatches to a specialized implementation based on the file type
Returns a tuple containing two vectors: the binned m/z axis and the summed intensities.
"""
-function get_total_spectrum(msi_data::MSIData; num_bins::Int=2000)
+function get_total_spectrum(msi_data::MSIData; num_bins::Int=2000, mask_path::Union{String, Nothing}=nothing)::Tuple{Vector{Float64}, Vector{Float64}, Int}
+ local masked_indices::Union{Set{Int}, Nothing} = nothing
+ if mask_path !== nothing
+ mask_matrix = load_and_prepare_mask(mask_path, msi_data.image_dims)
+ masked_indices = get_masked_spectrum_indices(msi_data, mask_matrix)
+ println("Calculating total spectrum for masked region (mask from: $(mask_path))")
+ end
+
if msi_data.source isa ImzMLSource
- return get_total_spectrum_imzml(msi_data, num_bins=num_bins)
+ return get_total_spectrum_imzml(msi_data, num_bins=num_bins, masked_indices=masked_indices)
else # MzMLSource
- return get_total_spectrum_mzml(msi_data, num_bins=num_bins)
+ return get_total_spectrum_mzml(msi_data, num_bins=num_bins, masked_indices=masked_indices)
end
end
@@ -715,14 +765,15 @@ This is effectively the total spectrum divided by the number of spectra.
Returns a tuple containing two vectors: the binned m/z axis and the averaged intensities.
"""
-function get_average_spectrum(msi_data::MSIData; num_bins::Int=2000)
- mz_bins, intensity_sum = get_total_spectrum(msi_data, num_bins=num_bins)
+function get_average_spectrum(msi_data::MSIData; num_bins::Int=2000, mask_path::Union{String, Nothing}=nothing)::Tuple{Vector{Float64}, Vector{Float64}}
+ mz_bins, intensity_sum, num_spectra_processed = get_total_spectrum(msi_data, num_bins=num_bins, mask_path=mask_path)
- if isempty(intensity_sum)
- return (mz_bins, intensity_sum)
+ if isempty(intensity_sum) || num_spectra_processed == 0
+ @warn "Cannot calculate average spectrum: no spectra were processed."
+ return (mz_bins, Float64[])
end
- num_spectra = length(msi_data.spectra_metadata)
- average_intensity = intensity_sum ./ num_spectra
+
+ average_intensity = intensity_sum ./ num_spectra_processed
return (mz_bins, average_intensity)
end
@@ -855,13 +906,16 @@ is significantly faster for bulk processing than reading spectra one by one.
- `data`: The `MSIData` object.
- `source`: The `ImzMLSource`.
"""
-function _iterate_uncompressed_fast(f::Function, data::MSIData, source::ImzMLSource)
+function _iterate_uncompressed_fast(f::Function, data::MSIData, source::ImzMLSource, indices_to_iterate::Union{AbstractVector{Int}, Nothing})
# Optimized path for uncompressed data using buffer reuse
max_points = maximum(meta -> meta.mz_asset.encoded_length, data.spectra_metadata)
mz_buffer = Vector{source.mz_format}(undef, max_points)
int_buffer = Vector{source.intensity_format}(undef, max_points)
- for i in 1:length(data.spectra_metadata)
+ # Determine which indices to iterate over
+ spectrum_indices = (indices_to_iterate === nothing) ? (1:length(data.spectra_metadata)) : indices_to_iterate
+
+ for i in spectrum_indices
meta = data.spectra_metadata[i]
nPoints = meta.mz_asset.encoded_length
@@ -904,10 +958,14 @@ and will be slower and allocate more memory than `_iterate_uncompressed_fast`.
- `data`: The `MSIData` object.
- `source`: The `ImzMLSource`.
"""
-function _iterate_compressed_fast(f::Function, data::MSIData, source::ImzMLSource)
+function _iterate_compressed_fast(f::Function, data::MSIData, source::ImzMLSource, indices_to_iterate::Union{AbstractVector{Int}, Nothing})
# Path for datasets containing at least one compressed spectrum.
# This path reads and decompresses each spectrum individually.
- for i in 1:length(data.spectra_metadata)
+
+ # Determine which indices to iterate over
+ spectrum_indices = (indices_to_iterate === nothing) ? (1:length(data.spectra_metadata)) : indices_to_iterate
+
+ for i in spectrum_indices
meta = data.spectra_metadata[i]
if meta.mz_asset.encoded_length == 0 && meta.int_asset.encoded_length == 0
@@ -937,7 +995,7 @@ This function acts as a dispatcher:
- If all spectra are uncompressed, it uses a highly optimized path that reuses pre-allocated
buffers to minimize memory allocations and overhead.
"""
-function _iterate_spectra_fast_impl(f::Function, data::MSIData, source::ImzMLSource)
+function _iterate_spectra_fast_impl(f::Function, data::MSIData, source::ImzMLSource, indices_to_iterate::Union{AbstractVector{Int}, Nothing})
if isempty(data.spectra_metadata)
return
end
@@ -947,9 +1005,9 @@ function _iterate_spectra_fast_impl(f::Function, data::MSIData, source::ImzMLSou
data.spectra_metadata)
if any_compressed
- _iterate_compressed_fast(f, data, source)
+ _iterate_compressed_fast(f, data, source, indices_to_iterate)
else
- _iterate_uncompressed_fast(f, data, source)
+ _iterate_uncompressed_fast(f, data, source, indices_to_iterate)
end
end
@@ -961,12 +1019,16 @@ through each spectrum, decodes the Base64 data on the fly, and calls the
provided function. It bypasses the `GetSpectrum` cache to avoid storing
all decoded spectra in memory.
"""
-function _iterate_spectra_fast_impl(f::Function, data::MSIData, source::MzMLSource)
+function _iterate_spectra_fast_impl(f::Function, data::MSIData, source::MzMLSource, indices_to_iterate::Union{AbstractVector{Int}, Nothing})
# This implementation is for mzML. It iterates through each spectrum,
# decodes it, and then calls the function `f`. It's less performant than
# the imzML version because we cannot pre-allocate buffers of a known size,
# but it is correct and still faster than using GetSpectrum due to bypassing the cache.
- for i in 1:length(data.spectra_metadata)
+
+ # Determine which indices to iterate over
+ spectrum_indices = (indices_to_iterate === nothing) ? (1:length(data.spectra_metadata)) : indices_to_iterate
+
+ for i in spectrum_indices
meta = data.spectra_metadata[i]
mz = read_binary_vector(source.file_handle, meta.mz_asset)
@@ -986,7 +1048,7 @@ It dispatches to a specialized implementation based on the data source type
- `f`: A function to call for each spectrum, with signature `f(index, mz, intensity)`.
- `data`: The `MSIData` object.
"""
-function _iterate_spectra_fast(f::Function, data::MSIData)
+function _iterate_spectra_fast(f::Function, data::MSIData, indices_to_iterate::Union{AbstractVector{Int}, Nothing}=nothing)
# Dispatch to the correct implementation based on the source type
- _iterate_spectra_fast_impl(f, data, data.source)
+ _iterate_spectra_fast_impl(f, data, data.source, indices_to_iterate)
end
\ No newline at end of file
diff --git a/src/MSI_src.jl b/src/MSI_src.jl
index 19c0cf6..72adaf3 100644
--- a/src/MSI_src.jl
+++ b/src/MSI_src.jl
@@ -14,7 +14,10 @@ export OpenMSIData,
get_average_spectrum,
LoadMzml,
precompute_analytics,
- process_spectrum
+ process_spectrum,
+ generate_colorbar_image,
+ process_image_pipeline,
+ load_and_prepare_mask
# Export the public Preprocessing API
export FeatureMatrix,
@@ -43,6 +46,7 @@ include("mzML.jl")
include("imzML.jl")
include("MzmlConverter.jl")
include("Preprocessing.jl")
+include("ImageProcessing.jl")
# --- Main Entry Point --- #
diff --git a/src/imzML.jl b/src/imzML.jl
index 6c3b48a..635eefa 100644
--- a/src/imzML.jl
+++ b/src/imzML.jl
@@ -744,10 +744,17 @@ This is a performant function that iterates through spectra once.
# Returns
- A `Matrix{Float64}` representing the intensity slice.
"""
-function get_mz_slice(data::MSIData, mass::Real, tolerance::Real)
+function get_mz_slice(data::MSIData, mass::Real, tolerance::Real; mask_path::Union{String, Nothing}=nothing)
width, height = data.image_dims
slice_matrix = zeros(Float64, height, width)
+ local masked_indices::Union{Set{Int}, Nothing} = nothing
+ if mask_path !== nothing
+ mask_matrix = load_and_prepare_mask(mask_path, (width, height))
+ masked_indices = get_masked_spectrum_indices(data, mask_matrix)
+ println("Applying mask from: $(mask_path), found $(length(masked_indices)) spectra in ROI.")
+ end
+
# INTELLIGENT LOADING: Ensure analytics are computed for filtering.
if data.spectrum_stats_df === nothing || !hasproperty(data.spectrum_stats_df, :MinMZ)
println("Per-spectrum metadata not found. Running one-time analytics computation...")
@@ -761,7 +768,9 @@ function get_mz_slice(data::MSIData, mass::Real, tolerance::Real)
# 1. Find all candidate spectra first for efficient filtering
candidate_indices = Set{Int}()
- for i in 1:length(data.spectra_metadata)
+ indices_to_check = masked_indices === nothing ? (1:length(data.spectra_metadata)) : masked_indices
+
+ for i in indices_to_check
spec_min_mz = stats_df.MinMZ[i]
spec_max_mz = stats_df.MaxMZ[i]
if target_max >= spec_min_mz && target_min <= spec_max_mz
@@ -769,20 +778,17 @@ function get_mz_slice(data::MSIData, mass::Real, tolerance::Real)
end
end
- println("Found $(length(candidate_indices)) candidate spectra (filtered from $(length(data.spectra_metadata)))")
+ println("Found $(length(candidate_indices)) candidate spectra (filtered from $(length(indices_to_check)) initial spectra)")
# 2. Iterate using the optimized, low-allocation iterator
results_count = 0
- _iterate_spectra_fast(data) do idx, mz_array, intensity_array
- # Process only the spectra that are candidates
- if idx in candidate_indices
- intensity = find_mass(mz_array, intensity_array, mass, tolerance)
- if intensity > 0.0
- meta = data.spectra_metadata[idx]
- if 1 <= meta.x <= width && 1 <= meta.y <= height
- slice_matrix[meta.y, meta.x] = intensity
- results_count += 1
- end
+ _iterate_spectra_fast(data, collect(candidate_indices)) do idx, mz_array, intensity_array
+ meta = data.spectra_metadata[idx]
+ intensity = find_mass(mz_array, intensity_array, mass, tolerance)
+ if intensity > 0.0
+ if 1 <= meta.x <= width && 1 <= meta.y <= height
+ slice_matrix[meta.y, meta.x] = intensity
+ results_count += 1
end
end
end
@@ -802,7 +808,7 @@ This is a highly performant function that iterates through the full dataset only
# Returns
- A `Dict{Real, Matrix{Float64}}` mapping each mass to its intensity slice matrix.
"""
-function get_multiple_mz_slices(data::MSIData, masses::AbstractVector{<:Real}, tolerance::Real)
+function get_multiple_mz_slices(data::MSIData, masses::AbstractVector{<:Real}, tolerance::Real; mask_path::Union{String, Nothing}=nothing)
width, height = data.image_dims
# 1. Initialize a dictionary to hold the output slice matrices
@@ -811,6 +817,13 @@ function get_multiple_mz_slices(data::MSIData, masses::AbstractVector{<:Real}, t
slice_dict[mass] = zeros(Float64, height, width)
end
+ local masked_indices::Union{Set{Int}, Nothing} = nothing
+ if mask_path !== nothing
+ mask_matrix = load_and_prepare_mask(mask_path, (width, height))
+ masked_indices = get_masked_spectrum_indices(data, mask_matrix)
+ println("Applying mask from: $(mask_path), found $(length(masked_indices)) spectra in ROI.")
+ end
+
# 2. Ensure analytics are computed for filtering.
if data.spectrum_stats_df === nothing || !hasproperty(data.spectrum_stats_df, :MinMZ)
println("Per-spectrum metadata not found. Running one-time analytics computation...")
@@ -820,12 +833,13 @@ function get_multiple_mz_slices(data::MSIData, masses::AbstractVector{<:Real}, t
println("Filtering candidate spectra for $(length(masses)) m/z values...")
stats_df = data.spectrum_stats_df
candidate_indices = Set{Int}()
+ indices_to_check = masked_indices === nothing ? (1:length(data.spectra_metadata)) : masked_indices
# 3. Find all spectra that could contain *any* of the requested masses.
for mass in masses
target_min = mass - tolerance
target_max = mass + tolerance
- for i in 1:length(data.spectra_metadata)
+ for i in indices_to_check
# If already a candidate, no need to check again
if i in candidate_indices
continue
@@ -841,20 +855,17 @@ function get_multiple_mz_slices(data::MSIData, masses::AbstractVector{<:Real}, t
println("Found $(length(candidate_indices)) total candidate spectra.")
# 4. Iterate through the data a single time using the optimized iterator.
- _iterate_spectra_fast(data) do idx, mz_array, intensity_array
- # Process only the spectra that are candidates
- if idx in candidate_indices
- meta = data.spectra_metadata[idx]
- # For this single spectrum, check all masses of interest
- for mass in masses
- # Check if this spectrum's range actually covers the current mass
- # This is a finer-grained check than the initial filtering
- if !isempty(mz_array) && (mass + tolerance) >= first(mz_array) && (mass - tolerance) <= last(mz_array)
- intensity = find_mass(mz_array, intensity_array, mass, tolerance)
- if intensity > 0.0
- if 1 <= meta.x <= width && 1 <= meta.y <= height
- slice_dict[mass][meta.y, meta.x] = intensity
- end
+ _iterate_spectra_fast(data, collect(candidate_indices)) do idx, mz_array, intensity_array
+ meta = data.spectra_metadata[idx]
+ # For this single spectrum, check all masses of interest
+ for mass in masses
+ # Check if this spectrum's range actually covers the current mass
+ # This is a finer-grained check than the initial filtering
+ if !isempty(mz_array) && (mass + tolerance) >= first(mz_array) && (mass - tolerance) <= last(mz_array)
+ intensity = find_mass(mz_array, intensity_array, mass, tolerance)
+ if intensity > 0.0
+ if 1 <= meta.x <= width && 1 <= meta.y <= height
+ slice_dict[mass][meta.y, meta.x] = intensity
end
end
end
@@ -878,30 +889,15 @@ Generates and saves a plot of a single image slice for a given m/z value.
This function closely imitates the logic of the original `GetSlice` but uses
the modern `MSIData` access patterns and robust peak finding.
"""
-function plot_slice(msi_data::MSIData, mass::Real, tolerance::Real, output_dir::String; stage_name="slice_mz_$(mass)", bins=256)
+function plot_slice(msi_data::MSIData, mass::Real, tolerance::Real, output_dir::String;
+ stage_name="slice_mz_$(mass)", bins=256, mask_path::Union{String, Nothing}=nothing)
- # 1. Create an empty image for the slice, with dimensions matching plotting expectations
- width, height = msi_data.image_dims
- slice_matrix = zeros(Float64, height, width)
-
- # 2. Iterate through each spectrum to build the slice
+ # 1. Generate the slice, applying the mask if provided.
println("Generating slice for m/z $mass...")
- _iterate_spectra_fast(msi_data) do spec_idx, mz_array, intensity_array
- meta = msi_data.spectra_metadata[spec_idx]
-
- # Find the peak intensity using the modern, robust find_mass
- intensity = find_mass(mz_array, intensity_array, mass, tolerance)
-
- if intensity > 0.0
- # Populate the matrix using (y, x) indexing
- if 1 <= meta.x <= width && 1 <= meta.y <= height
- slice_matrix[meta.y, meta.x] = intensity
- end
- end
- end
+ slice_matrix = get_mz_slice(msi_data, mass, tolerance, mask_path=mask_path)
println("Slice generation complete.")
- # 3. Plot the resulting slice matrix using CairoMakie
+ # 2. Plot the resulting slice matrix using CairoMakie
println("Plotting slice...")
fig = Figure(size = (600, 500))
@@ -925,7 +921,7 @@ function plot_slice(msi_data::MSIData, mass::Real, tolerance::Real, output_dir::
Colorbar(fig[1, 2], hm, label="Intensity")
colgap!(fig.layout, 5)
- # 4. Save the plot
+ # 3. Save the plot
mkpath(output_dir)
filename = "$(stage_name).png"
save_path = joinpath(output_dir, filename)
@@ -957,9 +953,9 @@ is used to exclude outliers before normalization.
# Returns
- A tuple `(low, high)` representing the calculated lower and upper intensity bounds.
"""
-function get_outlier_thres(img::AbstractMatrix{<:Real}, prob::Real=0.98)
+function get_outlier_thres(img::AbstractVector{<:Real}, prob::Real=0.98)
# DO NOT filter zeros. Use all pixel values like R does.
- int_values = vec(img)
+ int_values = img
low = minimum(int_values)
upp = maximum(int_values)
@@ -1001,6 +997,10 @@ function get_outlier_thres(img::AbstractMatrix{<:Real}, prob::Real=0.98)
return (low, actual_threshold)
end
+function get_outlier_thres(img::AbstractMatrix{<:Real}, prob::Real=0.98)
+ return get_outlier_thres(vec(img), prob)
+end
+
"""
set_pixel_depth(img, bounds, depth)
@@ -1063,11 +1063,26 @@ within these bounds.
# Returns
- A new image matrix with intensities quantized to the specified depth within the TrIQ bounds.
"""
-function TrIQ(pixMap::AbstractMatrix{<:Real}, depth::Integer, prob::Real=0.98)
- # Compute new dynamic range
- bounds = get_outlier_thres(pixMap, prob)
+function TrIQ(pixMap::AbstractMatrix{<:Real}, depth::Integer, prob::Real=0.98; mask_matrix::Union{BitMatrix, Nothing}=nothing)
+ local values_for_thres
+ if mask_matrix !== nothing
+ # Extract values only from the masked region for threshold calculation
+ values_for_thres = pixMap[mask_matrix]
+ # Filter out only zeros created by the mask.
+ filter!(x -> x != 0, values_for_thres)
+ else
+ values_for_thres = pixMap
+ end
+
+ # If the relevant values are empty or all zero, return a zero matrix
+ if isempty(values_for_thres) || all(iszero, values_for_thres)
+ bounds = (0.0, 0.0)
+ else
+ # Compute new dynamic range from the (potentially filtered) values
+ bounds = get_outlier_thres(values_for_thres, prob)
+ end
- # Set intensity dynamic range
+ # Set intensity dynamic range for the *entire* pixMap, using bounds derived from masked data
return set_pixel_depth(pixMap, bounds, depth)
end
@@ -1107,8 +1122,17 @@ Linearly scales the intensity values in a slice to a specified number of levels.
The output is an array of `UInt8` values.
This is a modernized version of `IntQuantCl`.
"""
-function quantize_intensity(slice::AbstractMatrix{<:Real}, levels::Integer=256)
- max_val = maximum(slice)
+function quantize_intensity(slice::AbstractMatrix{<:Real}, levels::Integer=256; mask_matrix::Union{BitMatrix, Nothing}=nothing)
+ local max_val
+ if mask_matrix !== nothing
+ masked_slice = slice[mask_matrix]
+ # Filter out zeros that might have been introduced by masking, if any
+ filter!(x -> x != 0, masked_slice)
+ max_val = isempty(masked_slice) ? 0.0 : maximum(masked_slice)
+ else
+ max_val = maximum(slice)
+ end
+
if max_val <= 0
return zeros(UInt8, size(slice))
end
@@ -1208,7 +1232,7 @@ function display_statistics(slices::AbstractArray{<:AbstractMatrix{<:Real}, 2},
stats_to_calc = Dict{String, Function}("Mean" => mean)
all_dfs = Dict{String, DataFrame}()
- for (stat_name, stat_func) in stats_to_to_calc
+ for (stat_name, stat_func) in stats_to_calc
# Pre-allocate matrix for better performance
stat_matrix = zeros(Float64, n_files, n_masses)
@@ -1497,4 +1521,81 @@ function process_columns_first!(A::AbstractMatrix)
end
# Viridis color palette (256 colors) - defined as constant
-const ViridisPalette = generate_palette(ColorSchemes.viridis)
\ No newline at end of file
+const ViridisPalette = generate_palette(ColorSchemes.viridis)
+
+function generate_colorbar_image(slice_data::AbstractMatrix, color_levels::Int, output_path::String; use_triq::Bool=false, triq_prob::Float64=0.98, mask_path::Union{String, Nothing}=nothing)
+ # 1. Determine bounds based on whether TrIQ is used and if a mask is applied
+ local data_for_bounds
+ if mask_path !== nothing
+ height, width = size(slice_data)
+ mask_matrix = load_and_prepare_mask(mask_path, (width, height))
+ # Extract only the non-zero values within the mask to accurately calculate the range
+ data_for_bounds = slice_data[mask_matrix]
+ filter!(x -> x > 0, data_for_bounds)
+ else
+ data_for_bounds = vec(slice_data)
+ end
+
+ # Handle cases where data_for_bounds might be empty after filtering
+ if isempty(data_for_bounds)
+ min_val, max_val = 0.0, 1.0 # Default range if no data in ROI
+ else
+ min_val, max_val = if use_triq
+ get_outlier_thres(data_for_bounds, triq_prob)
+ else
+ extrema(data_for_bounds)
+ end
+ end
+
+ # 2. Replicate the tick calculation logic from plot_slices
+ bins = color_levels
+ 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]
+
+ # 3. Create and save the colorbar image
+ fig = Figure(size=(150, 250))
+ Colorbar(fig[1, 1],
+ colormap=cgrad(:viridis, bins, categorical=true),
+ # limits=(min_val, max_val),
+ limits=(levels[1], levels[end]),
+ label=(scale == 1 ? "Intensity" : "Intensity ×10^$(round(Int, 3 * exponent))"),
+ ticks=(tick_positions, tick_labels),
+ labelsize=20,
+ ticklabelsize=16
+ )
+ save(output_path, fig)
+end
-