From 8fa92dc0e52f5f200e21475da6d362eb25856a64 Mon Sep 17 00:00:00 2001
From: Pixelguy14 <pixelguy14@gmail.com>
Date: Mon, 17 Nov 2025 13:31:37 -0600
Subject: [PATCH] improved preprocessing UI, added selective spectrum plot
 accordig to type (centroid, profile), potentially fixed a bug to retrieve
 folders on dropdown selectors in startup, improved preprocessing functions
 and test suite, testing a version control for preprocessing steps

---
 Manifest.toml                  |  97 +++-
 Project.toml                   |   7 +-
 app.jl                         |  33 +-
 app.jl.html                    | 392 ++++++++++----
 julia_imzML_visual.jl          |   6 +-
 mask.jl                        |  18 +
 mask.jl.html                   |   2 +-
 public/css/genieapp.css        |  33 +-
 src/BloomFilters.jl            |  64 ++-
 src/MSIData.jl                 |   6 +
 src/MSI_src.jl                 |   4 +
 src/Preprocessing.jl           | 927 ++++++++++++++++-----------------
 src/PreprocessingVersioning.jl | 446 ++++++++++++++++
 test/readme.md                 |   4 +-
 test/run_preprocessing.jl      | 735 ++++++++++++--------------
 15 files changed, 1758 insertions(+), 1016 deletions(-)
 create mode 100644 src/PreprocessingVersioning.jl

diff --git a/Manifest.toml b/Manifest.toml
index 62323b8..c6e5873 100644
--- a/Manifest.toml
+++ b/Manifest.toml
@@ -2,7 +2,12 @@
 
 julia_version = "1.11.7"
 manifest_format = "2.0"
-project_hash = "38f13d997585a9af185f1863ff626d584a7e791e"
+project_hash = "7843f141174fe17a3abc4f25c3819d3195459b5a"
+
+[[deps.ANSIColoredPrinters]]
+git-tree-sha1 = "574baf8110975760d391c710b6341da1afa48d8c"
+uuid = "a4c015fc-c6ff-483c-b24f-f7ea428134e9"
+version = "0.0.1"
 
 [[deps.ATK_jll]]
 deps = ["Artifacts", "Glib_jll", "JLLWrappers", "Libdl"]
@@ -159,6 +164,11 @@ git-tree-sha1 = "bca794632b8a9bbe159d56bf9e31c422671b35e0"
 uuid = "18cc8868-cbac-4acf-b575-c8ff214dc66f"
 version = "1.3.2"
 
+[[deps.Bessels]]
+git-tree-sha1 = "4435559dc39793d53a9e3d278e185e920b4619ef"
+uuid = "0e736298-9ec6-45e8-9647-e4fc86a2fe38"
+version = "0.2.8"
+
 [[deps.BitFlags]]
 git-tree-sha1 = "0691e34b3bb8be9307330f88d1a3c3f25466c24d"
 uuid = "d1d4a3ce-64b1-5f1a-9ba4-7e7e69966f35"
@@ -289,6 +299,12 @@ git-tree-sha1 = "980f01d6d3283b3dbdfd7ed89405f96b7256ad57"
 uuid = "da1fd8a2-8d9e-5ec2-8556-3022fb5608a2"
 version = "2.0.1"
 
+[[deps.CodecBase]]
+deps = ["TranscodingStreams"]
+git-tree-sha1 = "40956acdbef3d8c7cc38cba42b56034af8f8581a"
+uuid = "6c391c72-fb7b-5838-ba82-7cfb1bcfecbf"
+version = "0.3.4"
+
 [[deps.CodecZlib]]
 deps = ["TranscodingStreams", "Zlib_jll"]
 git-tree-sha1 = "962834c22b66e32aa10f7611c08c8ca4e20749a9"
@@ -397,6 +413,12 @@ weakdeps = ["IntervalSets", "LinearAlgebra", "StaticArrays"]
     ConstructionBaseLinearAlgebraExt = "LinearAlgebra"
     ConstructionBaseStaticArraysExt = "StaticArrays"
 
+[[deps.ContinuousWavelets]]
+deps = ["AbstractFFTs", "Documenter", "FFTW", "Interpolations", "LinearAlgebra", "SpecialFunctions", "Wavelets"]
+git-tree-sha1 = "6e883d34f7040ab2ec36c60c738ae53864819508"
+uuid = "96eb917e-2868-4417-9cb6-27e7ff17528f"
+version = "1.1.7"
+
 [[deps.Contour]]
 git-tree-sha1 = "439e35b0b36e2e5881738abc8857bd92ad6ff9a8"
 uuid = "d38c429a-6771-53c6-b99e-75d170b6e991"
@@ -434,6 +456,16 @@ git-tree-sha1 = "1a3f97f907e6dd8983b744d2642651bb162a3f7a"
 uuid = "dc8bdbbb-1ca9-579f-8c36-e416f6a65cce"
 version = "1.0.2"
 
+[[deps.DSP]]
+deps = ["Bessels", "FFTW", "IterTools", "LinearAlgebra", "Polynomials", "Random", "Reexport", "SpecialFunctions", "Statistics"]
+git-tree-sha1 = "5989debfc3b38f736e69724818210c67ffee4352"
+uuid = "717857b8-e6f2-59f4-9121-6e50c889abd2"
+version = "0.8.4"
+weakdeps = ["OffsetArrays"]
+
+    [deps.DSP.extensions]
+    OffsetArraysExt = "OffsetArrays"
+
 [[deps.DataAPI]]
 git-tree-sha1 = "abe83f3a2f1b857aac70ef8b269080af17764bbe"
 uuid = "9a962f9c-6df0-11e9-0e5d-c546b8b5ee8a"
@@ -516,6 +548,12 @@ git-tree-sha1 = "7442a5dfe1ebb773c29cc2962a8980f47221d76c"
 uuid = "ffbed154-4ef7-542d-bbb7-c09d3a79fcae"
 version = "0.9.5"
 
+[[deps.Documenter]]
+deps = ["ANSIColoredPrinters", "AbstractTrees", "Base64", "CodecZlib", "Dates", "DocStringExtensions", "Downloads", "Git", "IOCapture", "InteractiveUtils", "JSON", "Logging", "Markdown", "MarkdownAST", "Pkg", "PrecompileTools", "REPL", "RegistryInstances", "SHA", "TOML", "Test", "Unicode"]
+git-tree-sha1 = "352b9a04e74edd16429aec79f033620cf8e780d4"
+uuid = "e30172f5-a6a5-5a46-863b-614d45cd2de4"
+version = "1.15.0"
+
 [[deps.DotEnv]]
 deps = ["PrecompileTools"]
 git-tree-sha1 = "92e88cb68a5b10545234f46dfaeb2fa8a8a50c45"
@@ -651,12 +689,6 @@ git-tree-sha1 = "05882d6995ae5c12bb5f36dd2ed3f61c98cbb172"
 uuid = "53c48c17-4a7d-5ca2-90c5-79b7896eea93"
 version = "0.8.5"
 
-[[deps.FlameGraphs]]
-deps = ["AbstractTrees", "Colors", "FileIO", "FixedPointNumbers", "IndirectArrays", "LeftChildRightSiblingTrees", "Profile"]
-git-tree-sha1 = "0166baf81babb91cf78bfcc771d8e87c43d568df"
-uuid = "08572546-2f56-4bcf-ba4e-bab62c3a3f89"
-version = "1.1.0"
-
 [[deps.Fontconfig_jll]]
 deps = ["Artifacts", "Bzip2_jll", "Expat_jll", "FreeType2_jll", "JLLWrappers", "Libdl", "Libuuid_jll", "Zlib_jll"]
 git-tree-sha1 = "f85dac9a96a01087df6e3a749840015a0ca3817d"
@@ -804,6 +836,24 @@ git-tree-sha1 = "6570366d757b50fabae9f4315ad74d2e40c0560a"
 uuid = "59f7168a-df46-5410-90c8-f2779963d0ec"
 version = "5.2.3+0"
 
+[[deps.Git]]
+deps = ["Git_LFS_jll", "Git_jll", "JLLWrappers", "OpenSSH_jll"]
+git-tree-sha1 = "824a1890086880696fc908fe12a17bcf61738bd8"
+uuid = "d7ba0133-e1db-5d97-8f8c-041e4b3a1eb2"
+version = "1.5.0"
+
+[[deps.Git_LFS_jll]]
+deps = ["Artifacts", "JLLWrappers", "Libdl"]
+git-tree-sha1 = "bb8471f313ed941f299aa53d32a94ab3bee08844"
+uuid = "020c3dae-16b3-5ae5-87b3-4cb189e250b2"
+version = "3.7.0+0"
+
+[[deps.Git_jll]]
+deps = ["Artifacts", "Expat_jll", "JLLWrappers", "LibCURL_jll", "Libdl", "Libiconv_jll", "OpenSSL_jll", "PCRE2_jll", "Zlib_jll"]
+git-tree-sha1 = "b6a684587ebe896d9f68ae777f648205940f0f70"
+uuid = "f8c6e375-362e-5223-8a59-34ff63f689eb"
+version = "2.51.3+0"
+
 [[deps.Glib_jll]]
 deps = ["Artifacts", "GettextRuntime_jll", "JLLWrappers", "Libdl", "Libffi_jll", "Libiconv_jll", "Libmount_jll", "PCRE2_jll", "Zlib_jll"]
 git-tree-sha1 = "50c11ffab2a3d50192a228c313f05b5b5dc5acb2"
@@ -885,6 +935,12 @@ git-tree-sha1 = "68c173f4f449de5b438ee67ed0c9c748dc31a2ec"
 uuid = "34004b35-14d8-5ef3-9330-4cdb6864b03a"
 version = "0.3.28"
 
+[[deps.IOCapture]]
+deps = ["Logging", "Random"]
+git-tree-sha1 = "b6d6bfdd7ce25b0f9b2f6b3dd56b2673a66c8770"
+uuid = "b5f81e59-6552-4d32-b1f0-c071b021bf89"
+version = "0.2.5"
+
 [[deps.IfElse]]
 git-tree-sha1 = "debdd00ffef04665ccbb3e150747a77560e8fad1"
 uuid = "615f187c-cbe4-4ef1-ba3b-2fcf58d6d173"
@@ -1236,6 +1292,11 @@ git-tree-sha1 = "a9eaadb366f5493a5654e843864c13d8b107548c"
 uuid = "10f19ff3-798f-405d-979b-55457f8fc047"
 version = "0.1.17"
 
+[[deps.LazilyInitializedFields]]
+git-tree-sha1 = "0f2da712350b020bc3957f269c9caad516383ee0"
+uuid = "0e77f7df-68c5-4e49-93ce-4cd80f5598bf"
+version = "1.3.0"
+
 [[deps.LazyArtifacts]]
 deps = ["Artifacts", "Pkg"]
 uuid = "4af54fe1-eca0-43a8-85a7-787d91b784e3"
@@ -1436,6 +1497,12 @@ deps = ["Base64"]
 uuid = "d6f4376e-aef5-505a-96c1-9c027394607a"
 version = "1.11.0"
 
+[[deps.MarkdownAST]]
+deps = ["AbstractTrees", "Markdown"]
+git-tree-sha1 = "465a70f0fc7d443a00dcdc3267a497397b8a3899"
+uuid = "d0879d2d-cac2-40c8-9cee-1863dc0c7391"
+version = "0.1.2"
+
 [[deps.MathTeXEngine]]
 deps = ["AbstractTrees", "Automa", "DataStructures", "FreeTypeAbstraction", "GeometryBasics", "LaTeXStrings", "REPL", "RelocatableFolders", "UnicodeFun"]
 git-tree-sha1 = "a370fef694c109e1950836176ed0d5eabbb65479"
@@ -1785,10 +1852,6 @@ deps = ["Unicode"]
 uuid = "de0858da-6303-5e67-8744-51eddeeeb8d7"
 version = "1.11.0"
 
-[[deps.Profile]]
-uuid = "9abbd945-dff8-562f-b5e8-e1ebf5ef1b79"
-version = "1.11.0"
-
 [[deps.ProgressMeter]]
 deps = ["Distributed", "Printf"]
 git-tree-sha1 = "fbb92c6c56b34e1a2c4c36058f68f332bec840e7"
@@ -1872,6 +1935,12 @@ git-tree-sha1 = "4618ed0da7a251c7f92e869ae1a19c74a7d2a7f9"
 uuid = "dee08c22-ab7f-5625-9660-a9af2021b33f"
 version = "0.3.2"
 
+[[deps.RegistryInstances]]
+deps = ["LazilyInitializedFields", "Pkg", "TOML", "Tar"]
+git-tree-sha1 = "ffd19052caf598b8653b99404058fce14828be51"
+uuid = "2792f1a3-b283-48e8-9a74-f99dce5104f3"
+version = "0.1.0"
+
 [[deps.RelocatableFolders]]
 deps = ["SHA", "Scratch"]
 git-tree-sha1 = "ffdaf70d81cf6ff22c2b6e733c900c3321cab864"
@@ -2373,6 +2442,12 @@ git-tree-sha1 = "d1d9a935a26c475ebffd54e9c7ad11627c43ea85"
 uuid = "3d5dd08c-fd9d-11e8-17fa-ed2836048c2f"
 version = "0.21.72"
 
+[[deps.Wavelets]]
+deps = ["DSP", "FFTW", "LinearAlgebra", "Reexport", "SpecialFunctions", "Statistics"]
+git-tree-sha1 = "d0ec97a100abbe47a5e9a02361841da49cce6029"
+uuid = "29a6e085-ba6d-5f35-a997-948ac2efa89a"
+version = "0.10.1"
+
 [[deps.Wayland_jll]]
 deps = ["Artifacts", "EpollShim_jll", "Expat_jll", "JLLWrappers", "Libdl", "Libffi_jll"]
 git-tree-sha1 = "96478df35bbc2f3e1e791bc7a3d0eeee559e60e9"
diff --git a/Project.toml b/Project.toml
index 242e031..d21383f 100644
--- a/Project.toml
+++ b/Project.toml
@@ -7,12 +7,13 @@ Accessors = "7d9f7c33-5ae7-4f3b-8dc6-eff91059b697"
 Base64 = "2a0f44e3-6c83-55bd-87e4-b1978d98bd5f"
 CSV = "336ed68f-0bac-5ca0-87d4-7b16caf5d00b"
 CairoMakie = "13f3f980-e62b-5c42-98c6-ff1f3baf88f0"
+CodecBase = "6c391c72-fb7b-5838-ba82-7cfb1bcfecbf"
 ColorSchemes = "35d6a980-a343-548e-a6ea-1d62b119f2f4"
 Colors = "5ae59095-9a9b-59fe-a467-6f913c188581"
+ContinuousWavelets = "96eb917e-2868-4417-9cb6-27e7ff17528f"
 DataFrames = "a93c6f00-e57d-5684-b7b6-d8193f3e46c0"
 Dates = "ade2ca70-3891-5945-98fb-dc099432e06a"
 FileIO = "5789e2e9-d7fb-5bc7-8068-2c6fae9b9549"
-FlameGraphs = "08572546-2f56-4bcf-ba4e-bab62c3a3f89"
 GLMakie = "e9467ef8-e4e7-5192-8a1a-b1aee30e663a"
 Genie = "c43c736e-a2d1-11e8-161f-af95117fbd1e"
 GenieFramework = "a59fdf5c-6bf0-4f5d-949c-a137c9e2f353"
@@ -25,6 +26,8 @@ ImageFiltering = "6a3955dd-da59-5b1f-98d4-e7296123deb5"
 ImageMorphology = "787d08f9-d448-5407-9aad-5290dd7ab264"
 ImageSegmentation = "80713f31-8817-5129-9cf8-209ff8fb23e1"
 Images = "916415d5-f1e6-5110-898d-aaa5f9f070e0"
+Interpolations = "a98d9a8b-a2ab-59e6-89dd-64a1c18fca59"
+JLD2 = "033835bb-8acc-5ee8-8aae-3f567f8a3819"
 JSON = "682c06a0-de6a-54ab-a142-c8b1cf79cde6"
 Libz = "2ec943e9-cfe8-584d-b93d-64dcb6d567b7"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
@@ -47,7 +50,7 @@ Accessors = "0.1"
 Base64 = "1.11"
 CSV = "0.10"
 CairoMakie = "0.13"
-ColorSchemes = "3.31"
+ColorSchemes = "3.30"
 Colors = "0.12"
 DataFrames = "1.7"
 Dates = "1.11"
diff --git a/app.jl b/app.jl
index 5e4a5ef..8dd47c6 100644
--- a/app.jl
+++ b/app.jl
@@ -213,6 +213,7 @@ end
 
     # == Batch Processing & Registry Variables ==
     @private registry_init_done = false
+    @in refetch_folders = false
     @in selected_files = String[]
     @in available_folders = String[]
     @in image_available_folders = String[]
@@ -1919,6 +1920,37 @@ end
         end
     end
 
+    @onbutton refetch_folders begin
+        # Re-load registry and update folder lists
+        registry = load_registry(registry_path)
+        all_folders = sort(collect(keys(registry)), lt=natural)
+        img_folders = filter(folder -> get(get(registry, folder, Dict()), "is_imzML", false), all_folders)
+        
+        available_folders = deepcopy(all_folders)
+        image_available_folders = deepcopy(img_folders)
+
+        # For q-selects using image_available_folders
+        if !isempty(image_available_folders)
+            first_img_folder = first(image_available_folders)
+            if isempty(selected_folder_main)
+                selected_folder_main = first_img_folder
+            end
+            if isempty(selected_folder_compare_left)
+                selected_folder_compare_left = first_img_folder
+            end
+            if isempty(selected_folder_compare_right)
+                selected_folder_compare_right = first_img_folder
+            end
+        end
+
+        # For q-selects using available_folders
+        if !isempty(available_folders)
+            if isempty(selected_folder_metadata)
+                selected_folder_metadata = first(available_folders)
+            end
+        end
+    end
+
     @mounted watchplots()
 
     @onchange isready @time begin
@@ -1967,7 +1999,6 @@ end
                 image_available_folders = deepcopy(img_folders)
                 
                 println("UI lists updated. All: $(length(available_folders)), Images: $(length(image_available_folders))")
-
             catch e
                 @warn "Registry synchronization failed: $e"
                 available_folders = []
diff --git a/app.jl.html b/app.jl.html
index 1aea080..47a3673 100644
--- a/app.jl.html
+++ b/app.jl.html
@@ -19,7 +19,7 @@
           <div class="text-h6">imzML & mzML Data Pre-Treatment</div>
           <div class="row items-center">
             <q-input standout="custom-standout" class="q-ma-sm cursor-pointer col" v-model="full_route" readonly
-              :label="batch_file_count > 0 ? batch_file_count + ' file(s) in batch' : 'Select an imzML file'"
+              :label="batch_file_count > 0 ? batch_file_count + ' file(s) in batch' : 'Select an imzMl / mzML file'"
               v-on:click="btnSearch=true">
               <template v-slot:append>
                 <q-icon name="search" v-on:click="btnSearch=true" class="cursor-pointer" />
@@ -29,6 +29,22 @@
             <q-btn class="q-ma-sm" icon="clear" v-on:click="clear_batch_btn=true" :disable="batch_file_count === 0"
               label="Clear"></q-btn>
           </div>
+          <q-list bordered separator v-if="selected_files.length > 0">
+            <q-item v-for="(file, index) in selected_files" :key="index">
+              <q-item-section>
+                {{ file }}
+              </q-item-section>
+              <q-item-section side>
+                <q-btn flat round icon="delete" size="sm" v-on:click="selected_files.splice(index, 1)"></q-btn>
+              </q-item-section>
+            </q-item>
+          </q-list>
+          <q-input standout="custom-standout" class="q-ma-sm cursor-pointer col" v-model="full_route_cal" readonly
+            label="Select a calibration imzMl / mzML file" v-on:click="btnSearchCal=true">
+            <template v-slot:append>
+              <q-icon name="search" v-on:click="btnSearchCal=true" class="cursor-pointer" />
+            </template>
+          </q-input>-
           <br>
           <br>
           <q-tabs v-model="pre_tab" dense class="text-grey" indicator-color="primary" align="justify">
@@ -40,10 +56,176 @@
             <q-tab name="peak" label="Peak Detection"></q-tab>
             <q-tab name="binning" label="Binning"></q-tab>
           </q-tabs>
-          <q-separator></q-separator>
+          <q-separator />
           <q-tab-panels v-model="pre_tab" animated>
-            <q-tab-panel name="stabilization">
-              <p>HI!</p>
+            <q-tab-panel name="stabilization" class="q-pa-md">
+              <q-card class="q-mb-md">
+                <q-card-section>
+                  <div class="text-h6">Method</div>
+                </q-card-section>
+                <q-card-section>
+                  <q-radio v-model="stabilization_method" val="sqrt" label="SQRT" /><br>
+                  <q-radio v-model="stabilization_method" val="log" label="LOG" /><br>
+                  <q-radio v-model="stabilization_method" val="log2" label="LOG 2" /><br>
+                  <q-radio v-model="stabilization_method" val="log10" label="LOG 10" /><br>
+                </q-card-section>
+              </q-card>
+              <div class="row justify-end">
+                <q-btn :loading="progressPrep" class="q-ma-sm btn-style" icon="play_arrow" v-on:click="acceptStab=true"
+                  padding="lg" label="Accept" />
+                <q-btn :loading="progressPrep" class="q-ma-sm btn-style" icon="clear" v-on:click="undoPrep=true"
+                  padding="lg" label="Undo" />
+              </div>
+            </q-tab-panel>
+            <q-tab-panel name="smoothing">
+              <q-card class="q-mb-md">
+                <q-card-section>
+                  <div class="text-h6">Method</div>
+                </q-card-section>
+                <q-card-section>
+                  <q-radio v-model="smoothing_method" val="sg" label="Savitzky-Golay" /><br>
+                  <q-radio v-model="smoothing_method" val="ma" label="Moving Average" /><br>
+                </q-card-section>
+              </q-card>
+              <q-card>
+                <q-card-section>
+                  <div class="text-h6">Parameters</div>
+                </q-card-section>
+                <q-card-section>
+                  <q-input standout="custom-standout" label="Half Window Size" type="number"
+                    v-model.number="smoothing_window"></q-input>
+                </q-card-section>
+              </q-card>
+              <div class="row justify-end q-mt-md">
+                <q-btn :loading="progressPrep" class="q-ma-sm btn-style" icon="play_arrow" v-on:click="acceptSmoo=true"
+                  padding="lg" label="Accept" />
+                <q-btn :loading="progressPrep" class="q-ma-sm btn-style" icon="clear" v-on:click="undoPrep=true"
+                  padding="lg" label="Undo" />
+              </div>
+            </q-tab-panel>
+            <q-tab-panel name="baseline">
+              <q-card class="q-mb-md">
+                <q-card-section>
+                  <div class="text-h6">Method</div>
+                </q-card-section>
+                <q-card-section>
+                  <q-radio v-model="baseline_method" val="snip" label="SNIP" /><br>
+                  <q-radio v-model="baseline_method" val="tophat" label="TOP HAT" /><br>
+                  <q-radio v-model="baseline_method" val="convex_hull" label="CONVEX HULL" /><br>
+                  <q-radio v-model="baseline_method" val="median" label="MEDIAN" /><br>
+                </q-card-section>
+              </q-card>
+              <q-card>
+                <q-card-section>
+                  <div class="text-h6">Parameters</div>
+                </q-card-section>
+                <q-card-section>
+                  <q-input standout="custom-standout" label="Iterations" type="number"
+                    v-model.number="baseline_iterations"></q-input>
+                </q-card-section>
+              </q-card>
+              <div class="row justify-end q-mt-md">
+                <q-btn :loading="progressPrep" class="q-ma-sm btn-style" icon="play_arrow" v-on:click="acceptBase=true"
+                  padding="lg" label="Accept" />
+                <q-btn :loading="progressPrep" class="q-ma-sm btn-style" icon="clear" v-on:click="undoPrep=true"
+                  padding="lg" label="Undo" />
+              </div>
+            </q-tab-panel>
+            <q-tab-panel name="calibration">
+              <q-card class="q-mb-md">
+                <q-card-section>
+                  <div class="text-h6">Method</div>
+                </q-card-section>
+                <q-card-section>
+                  <q-radio v-model="calibration_method" val="tic" label="TIC" /><br>
+                  <q-radio v-model="calibration_method" val="pqn" label="PQN" /><br>
+                  <q-radio v-model="calibration_method" val="median" label="MEDIAN" /><br>
+                </q-card-section>
+              </q-card>
+              <div class="row justify-end q-mt-md">
+                <q-btn :loading="progressPrep" class="q-ma-sm btn-style" icon="play_arrow" v-on:click="acceptCali=true"
+                  padding="lg" label="Accept" />
+                <q-btn :loading="progressPrep" class="q-ma-sm btn-style" icon="clear" v-on:click="undoPrep=true"
+                  padding="lg" label="Undo" />
+              </div>
+            </q-tab-panel>
+            <q-tab-panel name="warping">
+              <q-card class="q-mb-md">
+                <q-card-section>
+                  <div class="text-h6">Method</div>
+                </q-card-section>
+                <q-card-section>
+                  <q-radio v-model="warping_method" val="snip" label="SNIP" /><br>
+                  <q-radio v-model="warping_method" val="tophat" label="TOP HAT" /><br>
+                  <q-radio v-model="warping_method" val="convex_hull" label="CONVEX HULL" /><br>
+                  <q-radio v-model="warping_method" val="median" label="MEDIAN" /><br>
+                </q-card-section>
+              </q-card>
+              <q-card>
+                <q-card-section>
+                  <div class="text-h6">Parameters</div>
+                </q-card-section>
+                <q-card-section>
+                  <q-input standout="custom-standout" label="Iterations" type="number"
+                    v-model.number="warping_iterations"></q-input>
+                </q-card-section>
+              </q-card>
+              <div class="row justify-end q-mt-md">
+                <q-btn :loading="progressPrep" class="q-ma-sm btn-style" icon="play_arrow" v-on:click="acceptWarp=true"
+                  padding="lg" label="Accept" />
+                <q-btn :loading="progressPrep" class="q-ma-sm btn-style" icon="clear" v-on:click="undoPrep=true"
+                  padding="lg" label="Undo" />
+              </div>
+            </q-tab-panel>
+            <q-tab-panel name="peak">
+              <p>Detect if profile or centroid to determine which elements to show.</p>
+              <q-card class="q-mb-md">
+                <q-card-section>
+                  <div class="text-h6">Method</div>
+                </q-card-section>
+                <q-card-section>
+                  <q-radio v-model="peak_method" val="mad" label="MAD" /><br>
+                  <q-radio v-model="peak_method" val="super_smoother" label="SUPER SMOOTHER" /><br>
+                </q-card-section>
+              </q-card>
+              <q-card>
+                <q-card-section>
+                  <div class="text-h6">Parameters</div>
+                </q-card-section>
+                <q-card-section>
+                  <q-input standout="custom-standout" label="Signal to Noise" type="number"
+                    v-model.number="peak_snr"></q-input>
+                  <q-input standout="custom-standout" label="Half Window Size" type="number" class="q-mt-md"
+                    v-model.number="peak_window"></q-input>
+                  <q-input standout="custom-standout" label="Intensity Threshold" type="number" class="q-mt-md"
+                    v-model.number="peak_threshold"></q-input>
+                </q-card-section>
+              </q-card>
+              <div class="row justify-end q-mt-md">
+                <q-btn :loading="progressPrep" class="q-ma-sm btn-style" icon="play_arrow" v-on:click="acceptPeak=true"
+                  padding="lg" label="Accept" />
+                <q-btn :loading="progressPrep" class="q-ma-sm btn-style" icon="clear" v-on:click="undoPrep=true"
+                  padding="lg" label="Undo" />
+              </div>
+            </q-tab-panel>
+            <q-tab-panel name="binning">
+              <q-card>
+                <q-card-section>
+                  <div class="text-h6">Parameters</div>
+                </q-card-section>
+                <q-card-section>
+                  <q-input standout="custom-standout" label="Tolerance" type="number"
+                    v-model.number="binning_tolerance"></q-input>
+                  <q-input standout="custom-standout" label="Frequency Threshold" type="number" class="q-mt-md"
+                    v-model.number="binning_threshold"></q-input>
+                </q-card-section>
+              </q-card>
+              <div class="row justify-end q-mt-md">
+                <q-btn :loading="progressPrep" class="q-ma-sm btn-style" icon="play_arrow" v-on:click="acceptBinn=true"
+                  padding="lg" label="Accept" />
+                <q-btn :loading="progressPrep" class="q-ma-sm btn-style" icon="clear" v-on:click="undoPrep=true"
+                  padding="lg" label="Undo" />
+              </div>
             </q-tab-panel>
           </q-tab-panels>
         </q-tab-panel>
@@ -56,7 +238,7 @@
           </p>
           <div class="row items-center">
             <q-input standout="custom-standout" class="q-ma-sm cursor-pointer col" v-model="full_route" readonly
-              :label="batch_file_count > 0 ? batch_file_count + ' file(s) in batch' : 'Select an imzML file'"
+              :label="batch_file_count > 0 ? batch_file_count + ' file(s) in batch' : 'Select an imzMl / mzML file'"
               v-on:click="btnSearch=true">
               <template v-slot:append>
                 <q-icon name="search" v-on:click="btnSearch=true" class="cursor-pointer" />
@@ -171,10 +353,9 @@
             </q-btn>
             <q-btn icon="zoom_out_map" class="q-ma-sm on-right btn-style" v-on:click="compareBtn=true" padding="sm"
               label="Compare"></q-btn>
-            <!--<q-btn class="q-ma-sm btn-style" :disable="btnMetadataDisable"
-                          v-on:click="showMetadataBtn=true" label="Show Metadata"></q-btn>-->
             <q-btn class="q-ma-sm btn-style" icon="edit" label="Mask Editor" href="/mask"></q-btn>
-            <q-btn class="q-ma-sm btn-style" v-on:click="showMetadataBtn=true" label="Show Metadata"></q-btn>
+            <q-btn class="q-ma-sm btn-style" icon="dashboard" v-on:click="showMetadataBtn=true"
+              label="Show Metadata"></q-btn>
             <div class="q-pa-md row items-center" v-show="progress">
               <q-spinner color="primary" size="2em" class="q-mr-sm"></q-spinner>
               <div class="text-caption">{{ progress_message }}</div>
@@ -238,7 +419,6 @@
           <div class="text-h6">mzML to imzML Converter</div>
           <p>Select the .mzML file and the corresponding .txt synchronization file to convert them into an .imzML/.ibd
             pair.</p>
-
           <q-input standout="custom-standout" class="q-ma-sm cursor-pointer" v-model="mzml_full_route" readonly
             label="Select your .mzML file" v-on:click="btnSearchMzml=true">
             <template v-slot:append>
@@ -268,96 +448,118 @@
   <div class="row col-6">
     <!-- Right DIV -->
     <div id="intDivStyle-right" class="st-col col-12 col-sm st-module">
-      <st-tabs id="tabHeader-main" :ids="tabIDs" :labels="tabLabels" v-model="selectedTab" no-arrows></st-tabs>
-      <q-tab-panels v-model="selectedTab">
-        <q-tab-panel name="tab0">
-          <!-- Content for Tab 0 -->
-          <h6>Image visualizer</h6>
-          <div class="row items-center">
-            <q-select v-model="selected_folder_main" :options="image_available_folders" label="Select Dataset"
-              class="q-ma-sm" style="min-width: 200px;"></q-select>
-            <q-space></q-space>
-            <q-btn icon="arrow_back" class="q-my-sm btn-style" v-on:click="imgMinus=true"></q-btn>
-            <q-btn icon="arrow_forward" class="q-my-sm on-right btn-style" v-on:click="imgPlus=true"></q-btn>
-          </div>
-          <!-- Image manager -->
-          <div id="image-container-normal" class="row st-col col-12">
-            <div class="col-10 q-pa-none q-ma-none">
-              <plotly id="plotImg" :data="plotdataImg" :layout="plotlayoutImg" class="q-pa-none q-ma-none sync_data"
-                @click="data_click"></plotly>
+      <div v-if="left_tab === 'pre_treatment'">
+        <div class="text-h6 q-mb-md">Spectrum View</div>
+        <q-card class="q-mb-md">
+          <q-card-section>
+            <div class="text-subtitle1">Before Preprocessing</div>
+          </q-card-section>
+          <q-card-section>
+            <plotly id="plotSpectraBefore" :data="plotdata" :layout="plotlayout" class="q-pa-none q-ma-none"></plotly>
+          </q-card-section>
+        </q-card>
+        <q-card>
+          <q-card-section>
+            <div class="text-subtitle1">After Preprocessing</div>
+          </q-card-section>
+          <q-card-section>
+            <plotly id="plotSpectraAfter" :data="plotdata_after" :layout="plotlayout_after" class="q-pa-none q-ma-none">
+            </plotly>
+          </q-card-section>
+        </q-card>
+      </div>
+      <div v-else>
+        <st-tabs id="tabHeader-main" :ids="tabIDs" :labels="tabLabels" v-model="selectedTab" no-arrows></st-tabs>
+        <q-tab-panels v-model="selectedTab">
+          <q-tab-panel name="tab0">
+            <!-- Content for Tab 0 -->
+            <h6>Image visualizer</h6>
+            <div class="row items-center">
+              <q-select v-model="selected_folder_main" :options="image_available_folders" label="Select Dataset"
+                class="q-ma-sm" style="min-width: 200px;" v-on:focus="refetch_folders = true"></q-select>
+              <q-space></q-space>
+              <q-btn icon="arrow_back" class="q-my-sm btn-style" v-on:click="imgMinus=true"></q-btn>
+              <q-btn icon="arrow_forward" class="q-my-sm on-right btn-style" v-on:click="imgPlus=true"></q-btn>
             </div>
-            <div class="col-2 q-pa-none q-ma-none">
-              <q-img id="colorbar-normal" class="q-ma-none q-pa-none" :src="colorbar"></q-img>
+            <!-- Image manager -->
+            <div id="image-container-normal" class="row st-col col-12">
+              <div class="col-10 q-pa-none q-ma-none">
+                <plotly id="plotImg" :data="plotdataImg" :layout="plotlayoutImg" class="q-pa-none q-ma-none sync_data"
+                  @click="data_click"></plotly>
+              </div>
+              <div class="col-2 q-pa-none q-ma-none">
+                <q-img id="colorbar-normal" class="q-ma-none q-pa-none" :src="colorbar"></q-img>
+              </div>
             </div>
-          </div>
-          <p v-html="msgimg"></p>
-        </q-tab-panel>
+            <p v-html="msgimg"></p>
+          </q-tab-panel>
 
-        <q-tab-panel name="tab1">
-          <!-- Content for Tab 1 -->
-          <h6>TrIQ visualizer</h6>
-          <div class="row items-center">
-            <q-select v-model="selected_folder_main" :options="image_available_folders" label="Select Dataset"
-              class="q-ma-sm" style="min-width: 200px;"></q-select>
-            <q-space></q-space>
-            <q-btn icon="arrow_back" class="q-my-sm btn-style" v-on:click="imgMinusT=true"></q-btn>
-            <q-btn icon="arrow_forward" class="q-my-sm on-right btn-style" v-on:click="imgPlusT=true"></q-btn>
-          </div>
-          <!-- Triq Image manager -->
-          <div id="image-container-triq" class="row st-col col-12">
-            <div class="col-10 q-pa-none q-ma-none ">
-              <plotly id="plotImgT" :data="plotdataImgT" :layout="plotlayoutImgT" class="q-pa-none q-ma-none sync_data"
-                @click="data_click"></plotly>
+          <q-tab-panel name="tab1">
+            <!-- Content for Tab 1 -->
+            <h6>TrIQ visualizer</h6>
+            <div class="row items-center">
+              <q-select v-model="selected_folder_main" :options="image_available_folders" label="Select Dataset"
+                class="q-ma-sm" style="min-width: 200px;" v-on:focus="refetch_folders = true"></q-select>
+              <q-space></q-space>
+              <q-btn icon="arrow_back" class="q-my-sm btn-style" v-on:click="imgMinusT=true"></q-btn>
+              <q-btn icon="arrow_forward" class="q-my-sm on-right btn-style" v-on:click="imgPlusT=true"></q-btn>
             </div>
-            <div class="col-2 q-pa-none q-ma-none ">
-              <q-img id="colorbar-triq" class="q-ma-none q-pa-none" :src="colorbarT"></q-img>
+            <!-- Triq Image manager -->
+            <div id="image-container-triq" class="row st-col col-12">
+              <div class="col-10 q-pa-none q-ma-none ">
+                <plotly id="plotImgT" :data="plotdataImgT" :layout="plotlayoutImgT"
+                  class="q-pa-none q-ma-none sync_data" @click="data_click"></plotly>
+              </div>
+              <div class="col-2 q-pa-none q-ma-none ">
+                <q-img id="colorbar-triq" class="q-ma-none q-pa-none" :src="colorbarT"></q-img>
+              </div>
             </div>
-          </div>
-          <p v-html="msgtriq"></p>
-        </q-tab-panel>
+            <p v-html="msgtriq"></p>
+          </q-tab-panel>
 
-        <q-tab-panel name="tab2">
-          <div class="row items-center">
-            <q-select v-model="selected_folder_main" :options="available_folders" label="Select Dataset" class="q-ma-sm"
-              style="min-width: 200px;"></q-select>
-            <q-btn-dropdown class="q-ma-sm btn-style" :loading="progressSpectraPlot" :disable="btnSpectraDisable"
-              label="Generate Spectra" icon="play_arrow">
-              <template v-slot:loading>
-                <q-spinner-hourglass class="on-left" />
-                Loading plot
-              </template>
+          <q-tab-panel name="tab2">
+            <div class="row items-center">
+              <q-select v-model="selected_folder_main" :options="available_folders" label="Select Dataset"
+                class="q-ma-sm" style="min-width: 200px;" v-on:focus="refetch_folders = true"></q-select>
+              <q-btn-dropdown class="q-ma-sm btn-style" :loading="progressSpectraPlot" :disable="btnSpectraDisable"
+                label="Generate Spectra" icon="play_arrow">
+                <template v-slot:loading>
+                  <q-spinner-hourglass class="on-left" />
+                  Loading plot
+                </template>
 
-              <q-list>
-                <q-item clickable v-close-popup v-on:click="createMeanPlot=true">
-                  <q-item-section>
-                    <q-item-label>Mean spectrum plot</q-item-label>
-                  </q-item-section>
-                </q-item>
-                <q-item clickable v-close-popup v-on:click="createSumPlot=true">
-                  <q-item-section>
-                    <q-item-label>Sum Spectrum plot</q-item-label>
-                  </q-item-section>
-                </q-item>
-                <q-item clickable v-close-popup v-on:click="createXYPlot=true">
-                  <q-item-section>
-                    <q-item-label>Spectrum plot (X,Y)</q-item-label>
-                  </q-item-section>
-                </q-item>
-              </q-list>
-            </q-btn-dropdown>
-          </div>
-          <plotly id="plotSpectra" :data="plotdata" :layout="plotlayout" class="q-pa-none q-ma-none"></plotly>
-        </q-tab-panel>
+                <q-list>
+                  <q-item clickable v-close-popup v-on:click="createMeanPlot=true">
+                    <q-item-section>
+                      <q-item-label>Mean spectrum plot</q-item-label>
+                    </q-item-section>
+                  </q-item>
+                  <q-item clickable v-close-popup v-on:click="createSumPlot=true">
+                    <q-item-section>
+                      <q-item-label>Sum Spectrum plot</q-item-label>
+                    </q-item-section>
+                  </q-item>
+                  <q-item clickable v-close-popup v-on:click="createXYPlot=true">
+                    <q-item-section>
+                      <q-item-label>Spectrum plot (X,Y)</q-item-label>
+                    </q-item-section>
+                  </q-item>
+                </q-list>
+              </q-btn-dropdown>
+            </div>
+            <plotly id="plotSpectra" :data="plotdata" :layout="plotlayout" class="q-pa-none q-ma-none"></plotly>
+          </q-tab-panel>
 
-        <q-tab-panel name="tab3">
-          <!-- Content for Tab 3 -->
-          <plotly id="plotTopo" :data="plotdataC" :layout="plotlayoutC" class="q-pa-none q-ma-none"></plotly>
-        </q-tab-panel>
-        <q-tab-panel name="tab4">
-          <!-- Content for Tab 4 -->
-          <plotly id="plot3d" :data="plotdata3d" :layout="plotlayout3d" class="q-pa-none q-ma-none"></plotly>
-        </q-tab-panel>
-      </q-tab-panels>
+          <q-tab-panel name="tab3">
+            <!-- Content for Tab 3 -->
+            <plotly id="plotTopo" :data="plotdataC" :layout="plotlayoutC" class="q-pa-none q-ma-none"></plotly>
+          </q-tab-panel>
+          <q-tab-panel name="tab4">
+            <!-- Content for Tab 4 -->
+            <plotly id="plot3d" :data="plotdata3d" :layout="plotlayout3d" class="q-pa-none q-ma-none"></plotly>
+          </q-tab-panel>
+        </q-tab-panels>
+      </div>
     </div>
   </div>
 </div>
@@ -395,7 +597,8 @@
         <div class="col-6">
           <div class="row items-center">
             <q-select v-model="selected_folder_compare_left" :options="image_available_folders"
-              label="Select Left Dataset" class="q-ma-sm" style="min-width: 200px;"></q-select>
+              label="Select Left Dataset" class="q-ma-sm" style="min-width: 200px;"
+              v-on:focus="refetch_folders = true"></q-select>
             <q-space></q-space>
             <st-tabs id="tabHeaderCompareLeft" :ids="CompTabIDsLeft" :labels="CompTabLabelsLeft"
               v-model="CompSelectedTabLeft"></st-tabs>
@@ -464,7 +667,8 @@
         <div class="col-6">
           <div class="row items-center">
             <q-select v-model="selected_folder_compare_right" :options="image_available_folders"
-              label="Select Right Dataset" class="q-ma-sm" style="min-width: 200px;"></q-select>
+              label="Select Right Dataset" class="q-ma-sm" style="min-width: 200px;"
+              v-on:focus="refetch_folders = true"></q-select>
             <q-space></q-space>
             <st-tabs id="tabHeaderCompareRight" :ids="CompTabIDsRight" :labels="CompTabLabelsRight"
               v-model="CompSelectedTabRight"></st-tabs>
@@ -561,7 +765,7 @@
       <div class="text-h6">Dataset Summary</div>
       <q-space />
       <q-select v-model="selected_folder_metadata" :options="available_folders" label="Select Dataset" class="q-ma-sm"
-        style="min-width: 250px;" standout="custom-standout"></q-select>
+        style="min-width: 250px;" standout="custom-standout" v-on:focus="refetch_folders = true"></q-select>
     </q-card-section>
 
     <q-card-section class="q-pt-none">
diff --git a/julia_imzML_visual.jl b/julia_imzML_visual.jl
index b9b40c4..cedb8d9 100644
--- a/julia_imzML_visual.jl
+++ b/julia_imzML_visual.jl
@@ -666,10 +666,10 @@ function xySpectrumPlot(data::MSIData, xCoord::Int, yCoord::Int, imgWidth::Int,
     # Downsample for plotting performance
     mz_down, int_down = MSI_src.downsample_spectrum(mz, intensity)
 
-    trace = if spectrum_mode == MSI_src.PROFILE
-        PlotlyBase.scatter(x=mz_down, y=int_down, mode="lines", marker=attr(size=1, color="blue", opacity=0.5), name="Spectrum", hoverinfo="x", hovertemplate="<b>m/z</b>: %{x:.4f}<extra></extra>")
-    else
+    trace = if spectrum_mode == MSI_src.CENTROID
         PlotlyBase.stem(x=mz_down, y=int_down, marker=attr(size=1, color="blue", opacity=0.5), name="Spectrum", hoverinfo="x", hovertemplate="<b>m/z</b>: %{x:.4f}<extra></extra>")
+    else
+        PlotlyBase.scatter(x=mz_down, y=int_down, mode="lines", marker=attr(size=1, color="blue", opacity=0.5), name="Spectrum", hoverinfo="x", hovertemplate="<b>m/z</b>: %{x:.4f}<extra></extra>")
     end
 
     plotdata = [trace]
diff --git a/mask.jl b/mask.jl
index c426bf4..45365a3 100644
--- a/mask.jl
+++ b/mask.jl
@@ -151,6 +151,7 @@ end
     @out available_folders = String[]
     @out image_available_folders = String[]
     @private registry_init_done = false
+    @in refetch_folders = false
 
     @out imgInt = "" # Path to current slice
     @out current_msi = ""
@@ -711,6 +712,23 @@ end
         end
     end
 
+    @onbutton refetch_folders begin
+        # First, re-run the logic to populate the folder list
+        registry = load_registry(registry_path)
+        all_folders = sort(collect(keys(registry)), lt=natural)
+        img_folders = filter(folder -> get(get(registry, folder, Dict()), "is_imzML", false), all_folders)
+        
+        available_folders = deepcopy(all_folders)
+        image_available_folders = deepcopy(img_folders)
+
+        # Apply the logic to select the first item if none is selected
+        if !isempty(image_available_folders) && isempty(selected_folder_main)
+            # This ensures the list is up-to-date
+            image_available_folders = deepcopy(img_folders)
+            selected_folder_main = first(image_available_folders)
+        end
+    end
+
     @onchange isready begin
         if isready && !registry_init_done
             sleep(1.0) # Give frontend time to initialize
diff --git a/mask.jl.html b/mask.jl.html
index 60d7011..2ae26f3 100644
--- a/mask.jl.html
+++ b/mask.jl.html
@@ -15,7 +15,7 @@
             <div class="text-subtitle1 q-mt-md">Step 1: Select Slice</div>
             <div class="row items-center">
                 <q-select v-model="selected_folder_main" :options="image_available_folders" label="Select Dataset"
-                    class="q-ma-sm col" standout="custom-standout" :disable="is_editing_mask"></q-select>
+                    class="q-ma-sm col" standout="custom-standout" :disable="is_editing_mask"  v-on:focus="refetch_folders = true"></q-select>
             </div>
             <p class="text-center" v-html="msgimg"></p>
 
diff --git a/public/css/genieapp.css b/public/css/genieapp.css
index 98ddaf2..9aa8eed 100644
--- a/public/css/genieapp.css
+++ b/public/css/genieapp.css
@@ -45,9 +45,10 @@
     color: #009f90;
 }
 
-#intDivStyle-left {
+#intDivStyle-left, #intDivStyle-right {
     border-radius: 10px;
     padding: 10px;
+    background-color: #fdfdfd; /* Add a background color to distinguish from the outer div */
 }
 
 /* Tab styling to match your theme */
@@ -81,13 +82,13 @@
     font-size: 0.9rem;
 }
 
-#intDivStyle-left .q-tab-panels {
-    height: 700px; /* Set this to accommodate your tallest content */
+#intDivStyle-left .q-tab-panels, #intDivStyle-right .q-tab-panels {
+    /* Removed fixed height */
 }
 
-#intDivStyle-left .q-tab-panel {
-    height: 100%;
-    overflow-y: auto; /* Add scroll if content overflows */
+#intDivStyle-left .q-tab-panel, #intDivStyle-right .q-tab-panel {
+    height: auto; /* Let content define height */
+    overflow-y: hidden; /* Remove scrollbar */
 }
 
 .pixelated-plot svg image {
@@ -98,4 +99,22 @@
 
 *{
     font-family: 'Roboto', 'Lato', sans-serif;
-}
\ No newline at end of file
+}
+
+.q-card {
+    border-radius: 8px !important;
+    box-shadow: 0 2px 4px rgba(0,0,0,0.1) !important;
+    margin-bottom: 16px !important;
+}
+
+.q-card__section--dark {
+    background: #f7f7f7 !important;
+}
+
+.q-radio__label {
+    font-size: 1rem !important;
+}
+
+.q-option-group > div {
+    margin-bottom: 8px;
+}
diff --git a/src/BloomFilters.jl b/src/BloomFilters.jl
index 2d95ed2..e5a6790 100644
--- a/src/BloomFilters.jl
+++ b/src/BloomFilters.jl
@@ -37,9 +37,20 @@ mutable struct StreamingBloomFilter
 end
 
 """
-    BloomFilter(expected_elements::Int, false_positive_rate::Float64=0.01; seed::UInt64=0x12345678)
+    BloomFilter{T}(expected_elements::Int, false_positive_rate::Float64=0.01; kwargs...) -> Return type
 
 Creates a Bloom filter optimized for the expected number of elements and desired false positive rate.
+
+# Arguments
+
+- `expected_elements::Int`: Argument description
+- `false_positive_rate::Float64`: Argument description
+    (**Default**: `0.01`)
+
+# Keywords
+
+- `seed::Union{UInt32,UInt64}`: Keyword description
+    (**Default**: `0x12345678`)
 """
 function BloomFilter{T}(expected_elements::Int, false_positive_rate::Float64=0.01; seed::Union{UInt32,UInt64}=0x12345678) where T
     # Convert seed to UInt64 for consistency
@@ -56,9 +67,12 @@ end
 """
     optimal_bit_size(n::Int, p::Float64) -> Int
 
-Calculates the optimal number of bits for a Bloom filter given:
-- n: expected number of elements
-- p: desired false positive rate
+Calculates the optimal number of bits for a Bloom filter
+
+# Arguments
+
+- `n::Int`: Expected number of elements
+- `p::Float64`: Desired false positive rate
 """
 function optimal_bit_size(n::Int, p::Float64)::Int
     if p <= 0.0 || p >= 1.0
@@ -73,6 +87,11 @@ end
     optimal_hash_count(n::Int, m::Int) -> Int
 
 Calculates the optimal number of hash functions for a Bloom filter.
+    
+# Arguments
+
+- `n::Int`: Expected number of elements
+- `p::Float64`: Desired false positive rate
 """
 function optimal_hash_count(n::Int, m::Int)::Int
     if n <= 0 || m <= 0
@@ -109,6 +128,11 @@ end
     Base.push!(bf::BloomFilter{T}, item::T)
 
 Adds an element to the Bloom filter.
+
+# Arguments
+
+- `bf::BloomFilter{T}`: Argument description
+- `item::T`: Argument description
 """
 function Base.push!(bf::BloomFilter{T}, item::T) where T
     hashes = hash_functions(item, bf.hash_count, bf.size, bf.seed)
@@ -169,27 +193,39 @@ function false_positive_rate(bf::BloomFilter)::Float64
 end
 
 """
-    fill_ratio(bf::BloomFilter) -> Float64
+    fill_ratio(bf::BloomFilter)::Float64 return count(bf.bits) / length(bf.bits) end -> Return type
 
 Returns the fraction of bits that are set to 1.
+
+# Arguments
+
+- `bf::BloomFilter`: Argument description
 """
 function fill_ratio(bf::BloomFilter)::Float64
     return count(bf.bits) / length(bf.bits)
 end
 
 """
-    is_empty(bf::BloomFilter) -> Bool
+    is_empty(bf::BloomFilter)::Bool return bf.count == 0 end -> Return type
 
 Checks if the Bloom filter is empty (no elements added).
+
+# Arguments
+
+- `bf::BloomFilter`: Argument description
 """
 function is_empty(bf::BloomFilter)::Bool
     return bf.count == 0
 end
 
 """
-    reset!(bf::BloomFilter)
+    reset!(bf::BloomFilter) fill!(bf.bits, false) bf.count = 0 return bf end -> Return type
 
 Clears the Bloom filter, removing all elements.
+
+# Arguments
+
+- `bf::BloomFilter`: Argument description
 """
 function reset!(bf::BloomFilter)
     fill!(bf.bits, false)
@@ -198,6 +234,20 @@ function reset!(bf::BloomFilter)
 end
 
 # Specialized constructor for empty Bloom filters
+"""
+    BloomFilter{T}(; kwargs...) -> Return type
+
+Description of the function
+
+# Keywords
+
+- `size::Int`: Keyword description
+    (**Default**: `100`)
+- `hash_count::Int`: Keyword description
+    (**Default**: `3`)
+- `seed::Union{UInt32,UInt64}`: Keyword description
+    (**Default**: `0x12345678`)
+"""
 function BloomFilter{T}(;size::Int=100, hash_count::Int=3, seed::Union{UInt32,UInt64}=0x12345678) where T
     seed_uint64 = UInt64(seed)
     bits = falses(size)
diff --git a/src/MSIData.jl b/src/MSIData.jl
index ed47a8b..6b7c8d7 100644
--- a/src/MSIData.jl
+++ b/src/MSIData.jl
@@ -420,6 +420,12 @@ function read_binary_vector(data::MSIData, io::IO, asset::SpectrumAsset)
     return out_array
 end
 
+function read_binary_vector(data::MSIData, ts_handle::ThreadSafeFileHandle, asset::SpectrumAsset)
+    lock(ts_handle.lock) do
+        return read_binary_vector(data, ts_handle.handle, asset)
+    end
+end
+
 # Overload for different source types
 """
     read_spectrum_from_disk(source::ImzMLSource, meta::SpectrumMetadata)
diff --git a/src/MSI_src.jl b/src/MSI_src.jl
index bf5df59..d878bd2 100644
--- a/src/MSI_src.jl
+++ b/src/MSI_src.jl
@@ -28,8 +28,12 @@ export FeatureMatrix,
        snip_baseline,
        tic_normalize,
        pqn_normalize,
+       median_normalize,
        detect_peaks_profile,
+       detect_peaks_wavelet,
+       detect_peaks_centroid,
        align_peaks_lowess,
+       find_calibration_peaks,
        bin_peaks,
        plot_stage_spectrum,
        calculate_ppm_error, 
diff --git a/src/Preprocessing.jl b/src/Preprocessing.jl
index 62e1e30..a1a8cff 100644
--- a/src/Preprocessing.jl
+++ b/src/Preprocessing.jl
@@ -16,6 +16,8 @@ using SavitzkyGolay # For SavitzkyGolay filtering
 using Dates         # For now()
 using CSV           # For writing CSV files
 using DataFrames    # For creating dataframes
+using ContinuousWavelets # For CWT peak detection
+using Interpolations # For calibration
 
 # =============================================================================
 # Data Structures
@@ -24,19 +26,29 @@ using DataFrames    # For creating dataframes
 """
     FeatureMatrix
 
-Structure to hold the final feature matrix, including m/z bin boundaries and sample indices.
-
-# Fields
-- `matrix`: The numerical matrix where rows are samples and columns are features (m/z bins).
-- `mz_bins`: A vector of tuples `(low_mz, high_mz)` for each feature column.
-- `sample_ids`: A vector of indices corresponding to the original spectra.
+A struct to hold the final feature matrix generated from the preprocessing pipeline.
 """
 struct FeatureMatrix
-    matrix::Array{Float64,2}                 # samples × features
-    mz_bins::Vector{Tuple{Float64,Float64}}  # [(low, high), ...]
+    matrix::Array{Float64,2}
+    mz_bins::Vector{Tuple{Float64,Float64}}
     sample_ids::Vector{Int}
 end
 
+"""
+    QCParameters
+
+A struct to hold comprehensive quality control parameters for validating spectra.
+"""
+struct QCParameters
+    min_tic_threshold::Float64      # Minimum total ion current
+    max_tic_threshold::Float64      # Maximum TIC (saturation check)
+    snr_threshold::Float64          # Minimum signal-to-noise for a spectrum to be considered valid
+    peak_count_range::Tuple{Int,Int} # Acceptable peak count range (min, max)
+    spatial_consistency_threshold::Float64 # For MSI spatial coherence (stub)
+    calibration_accuracy_ppm::Float64 # Maximum allowed PPM error for calibration lock masses
+end
+
+
 # =============================================================================
 # 0) Quality Control (QC)
 # =============================================================================
@@ -52,7 +64,7 @@ qc_is_empty(mz::AbstractVector, intensity::AbstractVector)::Bool =
 """
     qc_is_regular(mz) -> Bool
 
-Checks that the m/z axis is monotonically non-decreasing, as expected in a profile spectrum.
+Checks that the m/z axis is monotonically non-decreasing.
 """
 function qc_is_regular(mz::AbstractVector)
     n = length(mz)
@@ -73,26 +85,33 @@ end
     transform_intensity(intensity; method=:sqrt) -> Vector
 
 Applies a variance-stabilizing transformation to the intensity vector.
-Supported methods: `:sqrt` (default) and `:log1p`.
 """
 function transform_intensity(intensity::AbstractVector{<:Real}; method::Symbol=:sqrt)
     if method === :sqrt
         return sqrt.(max.(zero(eltype(intensity)), intensity))
     elseif method === :log1p
         return log1p.(max.(zero(eltype(intensity)), intensity))
+    elseif method === :log
+        return log.(max.(eps(eltype(intensity)), intensity))
+    elseif method === :log2
+        return log2.(max.(eps(eltype(intensity)), intensity))
+    elseif method === :log10
+        return log10.(max.(eps(eltype(intensity)), intensity))
     else
         return collect(float.(intensity))
     end
 end
 
 """
-    smooth_spectrum(y; window=21, order=2) -> Vector
+    smooth_spectrum(y; window=9, order=2) -> Vector
 
-Applies a Savitzky–Golay filter if `SavitzkyGolay.jl` is available.
-Otherwise, falls back to a simple (non-phase-correct) moving average.
+Applies a Savitzky–Golay filter to smooth the intensity data.
 """
 function smooth_spectrum(y::AbstractVector{<:Real}; window::Int=9, order::Int=2)
     win = isodd(window) ? window : window + 1
+    if length(y) < win
+        return y # Cannot smooth if data is smaller than window
+    end
     res = SavitzkyGolay.savitzky_golay(collect(float.(y)), win, order)
     return res.y
 end
@@ -104,22 +123,20 @@ end
 """
     snip_baseline(y, iterations=100) -> Vector
 
-Estimates the baseline using a simple 1D SNIP (Statistics-sensitive Non-linear
-Iterative Peak-clipping) algorithm. `iterations` controls the aggressiveness.
+Estimates the baseline of a spectrum using the SNIP algorithm.
 """
-function snip_baseline(y::AbstractVector{<:Real}, iterations::Int=100)
+function snip_baseline(y::AbstractVector{<:Real}; iterations::Int=100)
     n = length(y)
-    b = collect(float.(y))  # work copy
+    b = collect(float.(y))
     buf = similar(b)
     for k in 1:iterations
         copyto!(buf, b)
         @inbounds for i in 2:n-1
             buf[i] = min(b[i], 0.5 * (b[i-1] + b[i+1]))
         end
-        # Handle endpoints
         buf[1] = min(b[1], b[2])
         buf[end] = min(b[end], b[end-1])
-        b, buf = buf, b # Swap buffers
+        b, buf = buf, b
     end
     return b
 end
@@ -131,7 +148,7 @@ end
 """
     tic_normalize(y) -> Vector
 
-Normalizes intensities to the Total Ion Current (TIC). If the sum is zero, returns a copy.
+Normalizes spectrum intensities to the Total Ion Current (TIC).
 """
 function tic_normalize(y::AbstractVector{<:Real})
     s = sum(y)
@@ -141,108 +158,154 @@ end
 """
     pqn_normalize(M) -> Matrix
 
-Performs Probabilistic Quotient Normalization on a matrix `M` where columns are spectra.
+Performs Probabilistic Quotient Normalization (PQN) on a matrix of spectra.
 """
 function pqn_normalize(M::AbstractMatrix{<:Real})
     M_float = collect(float.(M))
-    # Calculate reference spectrum (median across all spectra)
     ref = mapslices(median, M_float; dims=2)[:,1]
-    
-    # Calculate quotients for each spectrum relative to the reference
     Q = similar(M_float)
     @inbounds for j in axes(M_float, 2)
         Q[:, j] = M_float[:, j] ./ (ref .+ eps(eltype(M_float)))
     end
-    
-    # Find the median quotient for each spectrum (scaling factor)
-    s = [median( @view Q[:, j]) for j in axes(Q, 2)]
-    
-    # Normalize the original matrix
+    s = [median(@view Q[:, j]) for j in axes(Q, 2)]
     @inbounds for j in axes(M_float, 2)
         M_float[:, j] ./= (s[j] + eps(eltype(M_float)))
     end
     return M_float
 end
 
+"""
+    median_normalize(y) -> Vector
+
+Normalizes spectrum intensities by dividing by the median intensity.
+"""
+function median_normalize(y::AbstractVector{<:Real})
+    m = median(y)
+    return m <= 0 ? collect(float.(y)) : collect(float.(y)) ./ m
+end
+
 # =============================================================================
-# 4) Peak Detection (for Profile Data)
+# 4) Peak Detection
 # =============================================================================
 
 """
-    detect_peaks_profile(mz, y; half_window=10, snr_threshold=2.0)
+    detect_peaks_profile(mz, y; ...)
 
-Detects local maxima with a signal-to-noise threshold (using MAD for noise estimation).
-Assumes profile-mode data and a monotonic m/z axis.
+Enhanced peak detection for profile-mode spectra with advanced filtering.
 """
-function detect_peaks_profile(mz::AbstractVector{<:Real},
-                              y::AbstractVector{<:Real};
+function detect_peaks_profile(mz::AbstractVector{<:Real}, y::AbstractVector{<:Real};
                               half_window::Int=10,
-                              snr_threshold::Float64=2.0)
+                              snr_threshold::Float64=2.0,
+                              min_peak_width::Real=0.01,
+                              max_peak_width::Real=2.0,
+                              peak_shape_threshold::Float64=0.7, # Stub
+                              merge_peaks_tolerance::Float64=0.002,
+                              min_peak_prominence::Float64=0.1)
     n = length(y)
     n < 3 && return (Float64[], Float64[])
     
-    # Noise estimation using Median Absolute Deviation (robust to peaks)
     noise_level = mad(y, normalize=true) + eps(Float64)
-
-    # Smooth the spectrum to make peak detection more robust
     ys = smooth_spectrum(y; window=max(5, 2*half_window+1), order=2)
 
-    peak_idx = Int[]
+    peak_indices = Int[]
     @inbounds for i in 2:n-1
         left = max(1, i - half_window)
         right = min(n, i + half_window)
         
-        local_max = ys[i]
-        # A point is a peak if it's the maximum in its neighborhood and above the SNR threshold
-        if local_max >= maximum( @view ys[left:right]) && (local_max > snr_threshold * noise_level)
-            # Ensure we only record one point for flat-topped peaks
-            if isempty(peak_idx) || (i - last(peak_idx) > half_window)
-                push!(peak_idx, i)
-            end
+        # Prominence check
+        prominence = ys[i] - max(minimum(@view ys[left:i]), minimum(@view ys[i:right]))
+        
+        if ys[i] >= maximum(@view ys[left:right]) && 
+           (ys[i] > snr_threshold * noise_level) &&
+           (prominence > min_peak_prominence * ys[i])
+            push!(peak_indices, i)
         end
     end
+    
+    # (STUB) Peak shape validation would be applied here
+    # e.g., fitting a Gaussian and checking R^2 > peak_shape_threshold
 
-    # Return original intensities at peak locations
-    pk_mz  = [float(mz[i]) for i in peak_idx]
-    pk_int = [float(y[i]) for i in peak_idx]
+    # Merge close peaks
+    if !isempty(peak_indices) && merge_peaks_tolerance > 0
+        merged_indices = [peak_indices[1]]
+        for i in 2:length(peak_indices)
+            if (mz[peak_indices[i]] - mz[last(merged_indices)]) > merge_peaks_tolerance
+                push!(merged_indices, peak_indices[i])
+            elseif y[peak_indices[i]] > y[last(merged_indices)]
+                merged_indices[end] = peak_indices[i] # Replace with more intense peak
+            end
+        end
+        peak_indices = merged_indices
+    end
+
+    # (STUB) Peak width filtering would be applied here
+    # This would require calculating FWHM for each peak, which is computationally intensive.
+
+    pk_mz  = [float(mz[i]) for i in peak_indices]
+    pk_int = [float(y[i]) for i in peak_indices]
     
     return (pk_mz, pk_int)
 end
 
 """
-    detect_peaks_centroid(mz, y; intensity_threshold=0.0)
+    detect_peaks_wavelet(mz, intensity; ...)
 
-Filters centroided data based on a minimum intensity threshold.
+Detects peaks using Continuous Wavelet Transform (CWT).
 """
-function detect_peaks_centroid(mz::AbstractVector{<:Real},
-                               y::AbstractVector{<:Real};
-                               intensity_threshold::Float64=0.0)
+function detect_peaks_wavelet(mz::AbstractVector, intensity::AbstractVector; scales=1:10, snr_threshold=3.0)
+    n = length(intensity)
+    n < 10 && return (Float64[], Float64[])
+    cwt_res = ContinuousWavelets.cwt(intensity, ContinuousWavelets.morl)
     
+    peak_indices = Int[]
+    noise_level = mad(intensity, normalize=true) + eps(Float64)
+
+    for i in 2:n-1
+        if intensity[i] > intensity[i-1] && intensity[i] > intensity[i+1] && intensity[i] > snr_threshold * noise_level
+            if abs(cwt_res[i, end]) > 0
+                push!(peak_indices, i)
+            end
+        end
+    end
+    
+    return (mz[peak_indices], intensity[peak_indices])
+end
+
+"""
+    detect_peaks_centroid(mz, y; ...)
+
+Filters peaks in centroid-mode data.
+"""
+function detect_peaks_centroid(mz::AbstractVector{<:Real}, y::AbstractVector{<:Real}; intensity_threshold::Float64=0.0)
     keep_indices = findall(y .>= intensity_threshold)
-    
     return (mz[keep_indices], y[keep_indices])
 end
 
 # =============================================================================
-# 5) Peak Alignment
+# 5) Peak Alignment & Calibration
 # =============================================================================
 
 """
-    align_peaks_lowess(ref_mz, tgt_mz; tolerance=0.002) -> warp::Function
+    align_peaks_lowess(ref_mz, tgt_mz; ...)
 
-Generates a warping function `warp(x)` to map target m/z values to reference m/z values.
-Uses a lightweight LOWESS-like approach with linear interpolation.
+Enhanced peak alignment with PPM tolerance and other constraints.
 """
-function align_peaks_lowess(ref_mz::Vector{<:Real},
-                            tgt_mz::Vector{<:Real};
-                            tolerance::Float64=0.002)
-    # Efficiently match peaks between sorted lists
-    pairs = Tuple{Float64,Float64}[]  # (target_mz, reference_mz)
+function align_peaks_lowess(ref_mz::Vector{<:Real}, tgt_mz::Vector{<:Real};
+                            tolerance::Float64=0.002,
+                            tolerance_unit::Symbol=:mz,
+                            max_shift_ppm::Float64=50.0,
+                            min_matched_peaks::Int=5)
+    pairs = Tuple{Float64,Float64}[]
     i = 1; j = 1
     while i <= length(tgt_mz) && j <= length(ref_mz)
+        tol = (tolerance_unit == :ppm) ? (ref_mz[j] * tolerance / 1e6) : tolerance
         dt = tgt_mz[i] - ref_mz[j]
-        if abs(dt) <= tolerance
-            push!(pairs, (float(tgt_mz[i]), float(ref_mz[j])))
+        
+        if abs(dt) <= tol
+            # Max shift check
+            if abs(dt) * 1e6 / ref_mz[j] <= max_shift_ppm
+                push!(pairs, (float(tgt_mz[i]), float(ref_mz[j])))
+            end
             i += 1; j += 1
         elseif dt < 0
             i += 1
@@ -251,42 +314,68 @@ function align_peaks_lowess(ref_mz::Vector{<:Real},
         end
     end
 
-    if length(pairs) < 3
-        @warn "Too few matching peaks for alignment. Returning identity function."
+    if length(pairs) < min_matched_peaks
+        @warn "Too few matching peaks ($(length(pairs)) < $min_matched_peaks). Returning identity function."
         return x -> float.(x)
     end
 
     t = [p[1] for p in pairs]
     r = [p[2] for p in pairs]
+    
+    # (STUB) A robust regression model (e.g., RANSAC) would be better here.
+    itp = linear_interpolation(t, r, extrapolation_bc=Line())
+    return itp
+end
 
-    # Lightly smooth the mapping to reduce noise
-    t_s = smooth_spectrum(t; window=5, order=2)
-    r_s = smooth_spectrum(r; window=5, order=2)
+"""
+    find_calibration_peaks(mz, intensity, reference_masses; ...)
 
-    # Return a function that performs linear interpolation for warping
-    function warp(x::AbstractVector{<:Real})
-        out = similar(collect(float.(x)))
-        for (k, xv) in enumerate(x)
-            if xv <= t_s[1]
-                # Linear extrapolation at the start
-                m = (r_s[2]-r_s[1]) / (t_s[2]-t_s[1] + eps())
-                out[k] = r_s[1] + m*(xv - t_s[1])
-            elseif xv >= t_s[end]
-                # Linear extrapolation at the end
-                m = (r_s[end]-r_s[end-1]) / (t_s[end]-t_s[end-1] + eps())
-                out[k] = r_s[end-1] + m*(xv - t_s[end-1])
-            else
-                # Linear interpolation for points in the middle
-                lo = searchsortedlast(t_s, xv)
-                hi = lo + 1
-                α = (xv - t_s[lo]) / (t_s[hi] - t_s[lo] + eps())
-                out[k] = (1-α)*r_s[lo] + α*r_s[hi]
-            end
+Finds peaks that match a list of reference masses.
+"""
+function find_calibration_peaks(mz::AbstractVector, intensity::AbstractVector, reference_masses::AbstractVector; ppm_tolerance=20.0)
+    matched_peaks = Dict{Float64, Float64}()
+    detected_mz, _ = detect_peaks_profile(mz, intensity)
+    
+    for ref_mass in reference_masses
+        tol = ref_mass * ppm_tolerance / 1e6
+        candidates = findall(m -> abs(m - ref_mass) <= tol, detected_mz)
+        if !isempty(candidates)
+            closest_peak_idx = argmin(abs.(detected_mz[candidates] .- ref_mass))
+            matched_peaks[ref_mass] = detected_mz[candidates[closest_peak_idx]]
         end
-        return out
     end
+    return matched_peaks
+end
 
-    return warp
+"""
+    calibrate_spectra(spectra, internal_standards; ...)
+
+Calibrates spectra using internal standards.
+"""
+function calibrate_spectra(spectra::Vector, internal_standards::Vector; ppm_tolerance=20.0)
+    calibrated_spectra = similar(spectra)
+    for (i, spec) in enumerate(spectra)
+        mz, intensity = spec[1], spec[2]
+
+        matched_peaks = find_calibration_peaks(mz, intensity, internal_standards; ppm_tolerance=ppm_tolerance)
+        if length(matched_peaks) < 2
+            @warn "Spectrum $i: Not enough calibration peaks found. Skipping."
+            calibrated_spectra[i] = spec
+            continue
+        end
+        measured = sort(collect(values(matched_peaks)))
+        theoretical = sort(collect(keys(matched_peaks)))
+        itp = linear_interpolation(measured, theoretical, extrapolation_bc=Line())
+        
+        new_mz = itp(mz)
+
+        if length(spec) == 3
+            calibrated_spectra[i] = (new_mz, intensity, spec[3])
+        else
+            calibrated_spectra[i] = (new_mz, intensity)
+        end
+    end
+    return calibrated_spectra
 end
 
 # =============================================================================
@@ -294,92 +383,199 @@ end
 # =============================================================================
 
 """
-    _find_bin_index(x, bins) -> Int
+    bin_peaks(all_pk_mz, all_pk_int, tolerance; ...)
 
-Efficiently finds the index of the bin `(low, high)` that contains `x` using binary search.
-Returns 0 if not found.
-"""
-function _find_bin_index(x::Float64, bins::Vector{Tuple{Float64,Float64}})
-    lo, hi = 1, length(bins)
-    while lo <= hi
-        mid = (lo + hi) >>> 1
-        b = bins[mid]
-        if x < b[1]
-            hi = mid - 1
-        elseif x > b[2]
-            lo = mid + 1
-        else
-            return mid
-        end
-    end
-    return 0
-end
-
-"""
-    bin_peaks(all_pk_mz, all_pk_int, tolerance; frequency_threshold=0.25)
-
-Groups peaks from all spectra into consensus m/z bins and creates a feature matrix.
-Filters out features that do not appear in a minimum fraction of spectra.
+Enhanced peak binning with adaptive and PPM-based parameters.
 """
 function bin_peaks(all_pk_mz::Vector{<:AbstractVector{<:Real}},
                    all_pk_int::Vector{<:AbstractVector{<:Real}},
-                   tolerance::Float64; frequency_threshold::Float64=0.25)
+                   tolerance::Float64; 
+                   frequency_threshold::Float64=0.25,
+                   tolerance_unit::Symbol=:mz,
+                   adaptive_tolerance::Bool=false, # Stub
+                   min_peak_per_bin::Int=2,
+                   max_bin_width_ppm::Float64=100.0,
+                   intensity_weighted_centers::Bool=true)
     ns = length(all_pk_mz)
     ns == 0 && return (zeros(0,0), Tuple{Float64,Float64}[])
 
-    # 1) Collect all unique peak m/z values and sort them
-    flat_mz = Float64[]
-    for v in all_pk_mz
-        append!(flat_mz, float.(v))
-    end
-    sort!(flat_mz)
-    isempty(flat_mz) && return (zeros(ns, 0), Tuple{Float64,Float64}[])
+    # Collect all peaks with their intensities and original spectrum index
+    all_peaks = [(float(m), float(i), s_idx) for s_idx in 1:ns for (m, i) in zip(all_pk_mz[s_idx], all_pk_int[s_idx])]
+    sort!(all_peaks, by=p->p[1])
+    
+    isempty(all_peaks) && return (zeros(ns, 0), Tuple{Float64,Float64}[])
 
-    # 2) Create contiguous m/z bins based on tolerance
-    bins = Tuple{Float64,Float64}[]
-    cur_lo = flat_mz[1]
-    cur_hi = flat_mz[1]
-    for x in @view flat_mz[2:end]
-        if x - cur_hi <= tolerance
-            cur_hi = x # Extend the current bin
+    # Create bins
+    bins = []
+    current_bin_peaks = [all_peaks[1]]
+    for p in @view all_peaks[2:end]
+        bin_center = mean(first.(current_bin_peaks))
+        tol = (tolerance_unit == :ppm) ? (bin_center * tolerance / 1e6) : tolerance
+        
+        if p[1] - last(current_bin_peaks)[1] <= tol
+            push!(current_bin_peaks, p)
         else
-            push!(bins, (cur_lo, cur_hi)) # Finalize old bin
-            cur_lo = x; cur_hi = x      # Start a new one
-        end
-    end
-    push!(bins, (cur_lo, cur_hi))
-
-    # 3) Create the feature matrix (samples x features) using max intensity per bin
-    X = zeros(Float64, ns, length(bins))
-    for i in 1:ns
-        for (mzv, iv) in zip(all_pk_mz[i], all_pk_int[i])
-            bidx = _find_bin_index(float(mzv), bins)
-            if bidx > 0
-                X[i, bidx] = max(X[i, bidx], float(iv))
+            if length(current_bin_peaks) >= min_peak_per_bin
+                push!(bins, current_bin_peaks)
             end
+            current_bin_peaks = [p]
+        end
+    end
+    length(current_bin_peaks) >= min_peak_per_bin && push!(bins, current_bin_peaks)
+
+    # Filter bins by max width
+    filter!(b -> (last(b)[1] - first(b)[1]) * 1e6 / mean(p[1] for p in b) <= max_bin_width_ppm, bins)
+
+    # Create feature matrix
+    X = zeros(Float64, ns, length(bins))
+    final_bins_boundaries = Vector{Tuple{Float64,Float64}}(undef, length(bins))
+    
+    for (j, bin_peaks) in enumerate(bins)
+        # Calculate bin center
+        local bin_center
+        if intensity_weighted_centers
+            weights = [p[2] for p in bin_peaks]
+            bin_center = sum(p[1]*p[2] for p in bin_peaks) / sum(weights)
+        else
+            bin_center = mean(p[1] for p in bin_peaks)
+        end
+        
+        final_bins_boundaries[j] = (first(bin_peaks)[1], last(bin_peaks)[1])
+
+        for p in bin_peaks
+            s_idx = p[3]
+            X[s_idx, j] = max(X[s_idx, j], p[2])
         end
     end
 
-    # 4) Filter features by minimum frequency
+    # Filter by frequency
     if frequency_threshold > 0
         present_count = vec(sum(X .> 0, dims=1))
         min_count = ceil(Int, frequency_threshold * ns)
         keep_mask = findall(present_count .>= min_count)
         X = X[:, keep_mask]
-        bins = bins[keep_mask]
+        final_bins_boundaries = final_bins_boundaries[keep_mask]
     end
 
-    return (X, bins)
+    return (X, final_bins_boundaries)
 end
 
+# =============================================================================
+# 7) Spatial & Advanced Processing (Stubs & New Functions)
+# =============================================================================
+
+"""
+    find_ppm_error_by_region(msi_data, region_masks, reference_peaks) -> Dict
+
+Calculates PPM error statistics for different spatial regions.
+`region_masks` is a Dict mapping region names (e.g., :tumor) to BitMatrix masks.
+"""
+function find_ppm_error_by_region(msi_data::MSIData, region_masks::Dict, reference_peaks::Dict)
+    regional_reports = Dict{Symbol, NamedTuple}()
+    
+    width, height = msi_data.image_dims
+    
+    for (region_name, mask) in region_masks
+        mask_height, mask_width = size(mask)
+        if mask_width != width || mask_height != height
+            @warn "Mask dimensions ($(mask_width)x$(mask_height)) for region '$region_name' do not match image dimensions ($(width)x$(height)). Skipping."
+            continue
+        end
+
+        indices = [
+            i for i in 1:length(msi_data.spectra_metadata) 
+            if msi_data.spectra_metadata[i].x > 0 &&
+               msi_data.spectra_metadata[i].y > 0 &&
+               mask[msi_data.spectra_metadata[i].y, msi_data.spectra_metadata[i].x]
+        ]
+        
+        if isempty(indices) continue end
+        
+        # Call analyze_mass_accuracy with the specific indices for the region
+        regional_reports[region_name] = analyze_mass_accuracy(msi_data, reference_peaks; spectrum_indices=indices)
+    end
+    return regional_reports
+end
+
+"""
+    regional_calibration(msi_data, region_masks, reference_peaks)
+
+(STUB) Applies different calibration models to different spatial regions.
+"""
+function regional_calibration(msi_data::MSIData, region_masks::Dict, reference_peaks::Dict)
+    @warn "regional_calibration is a stub and not fully implemented."
+    # 1. For each region in region_masks:
+    # 2. Create a region-specific calibration model using `calibrate_spectra`.
+    # 3. Apply the model to all spectra within that region.
+    # 4. Return a new MSIData object or modified spectra vector.
+    return msi_data # Return unmodified for now
+end
+
+# =============================================================================
+# 8) Advanced Peak Quality & Adaptive Parameters
+# =============================================================================
+
+"""
+    calculate_peak_quality_metrics(peak_mz, peak_intensity, local_spectrum) -> Dict
+
+(STUB) Calculates advanced quality metrics for a single peak.
+"""
+function calculate_peak_quality_metrics(peak_mz, peak_intensity, local_spectrum)
+    # A full implementation would calculate:
+    # - Sharpness (e.g., ratio of height to FWHM)
+    # - Symmetry (e.g., ratio of left/right half-widths)
+    # - Signal-to-Noise (using local noise estimation)
+    # - Isolation score (how close are other peaks)
+    return Dict(:sharpness => 1.0, :symmetry => 1.0, :snr => 10.0)
+end
+
+
+
+
+"""
+    calculate_adaptive_bin_tolerance(ppm_error_distribution) -> Float64
+
+Calculates an appropriate binning tolerance based on observed mass accuracy.
+"""
+function calculate_adaptive_bin_tolerance(ppm_error_distribution::Vector{Float64})
+    if isempty(ppm_error_distribution)
+        return 20.0 # Default if no data
+    end
+    # A robust strategy: mean + 3 * std deviation to capture ~99.7% of peaks
+    return mean(ppm_error_distribution) + 3 * std(ppm_error_distribution)
+end
+
+
 # =============================================================================
 # 7) Plotting Helper
 # =============================================================================
 
 """
-    plot_stage_spectrum(mz, intensity; title, ...)
+    plot_stage_spectrum(mz, intensity; title, xlabel, ylabel) -> Figure
 
-Returns a `CairoMakie.Figure` for a single spectrum trace. The caller is responsible for saving.
+A helper function to generate a plot of a single spectrum using `CairoMakie`. This is
+useful for visualizing the output of different preprocessing steps.
+
+# Arguments
+- `mz::AbstractVector`: The m/z vector for the x-axis.
+- `intensity::AbstractVector`: The intensity vector for the y-axis.
+
+# Keyword Arguments
+- `title::AbstractString`: The title of the plot.
+- `xlabel::AbstractString`: The label for the x-axis. Defaults to "m/z".
+- `ylabel::AbstractString`: The label for the y-axis. Defaults to "Intensity".
+
+# Returns
+- `CairoMakie.Figure`: A `Figure` object containing the plot. The caller is responsible for displaying or saving it.
+
+# Example
+```julia
+using CairoMakie
+mz = 100:0.1:110;
+intensity = rand(length(mz));
+fig = plot_stage_spectrum(mz, intensity, title="My Spectrum");
+save("my_spectrum.png", fig);
+```
 """
 function plot_stage_spectrum(mz::AbstractVector, intensity::AbstractVector;
                              title::AbstractString, xlabel::AbstractString="m/z",
@@ -395,255 +591,57 @@ function plot_stage_spectrum(mz::AbstractVector, intensity::AbstractVector;
     return fig
 end
 
-# =============================================================================
-# 8) Pipeline Orchestrator
-# =============================================================================
-
-"""
-    run_preprocessing_pipeline(spectra; steps, params, on_stage)
-
-Executes a flexible preprocessing pipeline on a vector of spectra.
-
-# Arguments
-- `spectra`: A vector of `(mz, intensity)` tuples.
-- `steps`: A vector of symbols defining the pipeline order (e.g., `[:qc, :smooth, :baseline, :peaks, :bin]`).
-- `params`: A dictionary of parameters for each step.
-- `on_stage`: An optional callback function `on_stage(stage_symbol; idx, mz, intensity)` executed after each step for logging or visualization.
-
-# Returns
-- A `FeatureMatrix` if `:bin` is in the steps, otherwise the vector of processed spectra.
-"""
-function run_preprocessing_pipeline(spectra::Vector;
-                                    steps::Vector{Symbol},
-                                    params::Dict=Dict(),
-                                    on_stage::Function=(;kwargs...)->nothing)
-
-    processed = deepcopy(spectra)  # Don't mutate the original input
-    reference_peaks = nothing      # For alignment
-
-    _emit(stage::Symbol, idx::Int, mz, y) = on_stage(stage; idx=idx, mz=mz, intensity=y)
-
-    for step in steps
-        @info "Running step: $step"
-
-        if step === :qc
-            for (i, (mz, y)) in enumerate(processed)
-                (isempty(mz) || isempty(y)) && continue
-                _emit(:qc_raw, i, mz, y)
-                qc_is_empty(mz, y) && @warn "Spectrum at index $i is empty."
-                !qc_is_regular(mz) && @warn "m/z axis at index $i is not monotonic."
-            end
-
-        elseif step === :transform
-            meth = get(params, :transform_method, :sqrt)
-            for i in eachindex(processed)
-                mz, y = processed[i]
-                y_new = transform_intensity(y; method=meth)
-                processed[i] = (mz, y_new)
-                _emit(:transform, i, mz, y_new)
-            end
-
-        elseif step === :smooth
-            win = get(params, :sg_window, 21)
-            ord = get(params, :sg_order, 2)
-            for i in eachindex(processed)
-                mz, y = processed[i]
-                y_smooth = smooth_spectrum(y; window=win, order=ord)
-                processed[i] = (mz, y_smooth)
-                _emit(:smooth, i, mz, y_smooth)
-            end
-
-        elseif step === :baseline
-            iters = get(params, :snip_iterations, 100)
-            for i in eachindex(processed)
-                mz, y = processed[i]
-                baseline = snip_baseline(y, iters)
-                y_corrected = max.(0.0, y .- baseline)
-                processed[i] = (mz, y_corrected)
-                _emit(:baseline, i, mz, y_corrected)
-            end
-
-        elseif step === :normalize
-            mode = get(params, :normalize_method, :tic)
-            if mode === :tic
-                for i in eachindex(processed)
-                    mz, y = processed[i]
-                    y_norm = tic_normalize(y)
-                    processed[i] = (mz, y_norm)
-                    _emit(:normalize, i, mz, y_norm)
-                end
-            elseif mode === :pqn
-                # Note: PQN assumes spectra are on a common m/z grid.
-                matrix = hcat([float.(p[2]) for p in processed]...)
-                matrix_norm = pqn_normalize(matrix)
-                for i in eachindex(processed)
-                    mz, _ = processed[i]
-                    processed[i] = (mz, view(matrix_norm, :, i))
-                    _emit(:normalize, i, mz, processed[i][2])
-                end
-            end
-
-        elseif step === :peaks
-            peak_results = Vector{Tuple{Vector{Float64},Vector{Float64}}}(undef, length(processed))
-            hw = get(params, :peak_half_window, 10)
-            snr = get(params, :peak_snr, 2.0)
-            for (i, (mz, y)) in enumerate(processed)
-                pk_mz, pk_int = detect_peaks_profile(mz, y; half_window=hw, snr_threshold=snr)
-                peak_results[i] = (pk_mz, pk_int)
-                # Emit with original mz axis but maybe stem plot of peaks?
-                _emit(:peaks, i, pk_mz, pk_int)
-            end
-            processed = peak_results
-            # Set reference for alignment
-            reference_peaks = isempty(processed) ? nothing : processed[1][1]
-
-        elseif step === :align
-            reference_peaks === nothing && (@error "Alignment requires a :peaks step first."; continue)
-            tol = get(params, :align_tolerance, 0.002)
-            for i in 2:length(processed)
-                tgt_peaks, intens = processed[i]
-                warp_func = align_peaks_lowess(reference_peaks, tgt_peaks; tolerance=tol)
-                processed[i] = (warp_func(tgt_peaks), intens)
-                _emit(:align, i, processed[i][1], processed[i][2])
-            end
-
-        elseif step === :bin
-            all_pks  = [s[1] for s in processed]
-            all_ints = [s[2] for s in processed]
-            tol  = get(params, :bin_tolerance, 0.002)
-            freq = get(params, :bin_min_frequency, 0.25)
-            mat, mz_bins = bin_peaks(all_pks, all_ints, tol; frequency_threshold=freq)
-            return FeatureMatrix(mat, mz_bins, collect(1:length(processed)))
-        end
-    end
-
-    @warn "Pipeline finished without a :bin step; returning processed spectra."
-    return processed
-end
-
-"""
-    run_preprocessing_pipeline(msi_data::MSIData, indices::Vector{Int}; steps, params, on_stage)
-
-Executes a flexible preprocessing pipeline on a subset of spectra from an MSIData object,
-with mode-aware logic for centroid and profile data.
-"""
-function run_preprocessing_pipeline(msi_data::MSIData, indices::Vector{Int};
-                                    steps::Vector{Symbol},
-                                    params::Dict=Dict(),
-                                    on_stage::Function=(;kwargs...)->nothing)
-
-    # This version of the pipeline is mode-aware.
-    # It processes spectra directly from the MSIData object.
-
-    processed_spectra = Vector{Tuple}(undef, length(indices))
-
-    # First, load all spectra and apply initial steps that run on individual spectra
-    for (i, spec_idx) in enumerate(indices)
-        mz, intensity = GetSpectrum(msi_data, spec_idx)
-        mode = msi_data.spectra_metadata[spec_idx].mode
-
-        on_stage(:qc_raw; idx=spec_idx, mz=mz, intensity=intensity)
-
-        for step in steps
-            if step === :transform
-                meth = get(params, :transform_method, :sqrt)
-                intensity = transform_intensity(intensity; method=meth)
-                on_stage(:transform; idx=spec_idx, mz=mz, intensity=intensity)
-            elseif step === :smooth && mode == PROFILE
-                win = get(params, :sg_window, 21)
-                ord = get(params, :sg_order, 2)
-                intensity = smooth_spectrum(intensity; window=win, order=ord)
-                on_stage(:smooth; idx=spec_idx, mz=mz, intensity=intensity)
-            elseif step === :baseline && mode == PROFILE
-                iters = get(params, :snip_iterations, 100)
-                baseline = snip_baseline(intensity, iters)
-                intensity = max.(0.0, intensity .- baseline)
-                on_stage(:baseline; idx=spec_idx, mz=mz, intensity=intensity)
-            elseif step === :normalize
-                norm_mode = get(params, :normalize_method, :tic)
-                if norm_mode === :tic
-                    intensity = tic_normalize(intensity)
-                    on_stage(:normalize; idx=spec_idx, mz=mz, intensity=intensity)
-                end
-            end
-        end
-        processed_spectra[i] = (mz, intensity, mode) # Store mode for peak detection
-    end
-
-    # Now, handle steps that require all spectra (like PQN) or are the final steps
-    final_result = nothing
-    for step in steps
-        if step === :normalize && get(params, :normalize_method, :tic) === :pqn
-            # Note: PQN assumes spectra are on a common m/z grid.
-            matrix = hcat([float.(p[2]) for p in processed_spectra]...)
-            matrix_norm = pqn_normalize(matrix)
-            for i in eachindex(processed_spectra)
-                mz, _, mode = processed_spectra[i]
-                processed_spectra[i] = (mz, view(matrix_norm, :, i), mode)
-                on_stage(:normalize; idx=indices[i], mz=mz, intensity=processed_spectra[i][2])
-            end
-        elseif step === :peaks
-            peak_results = Vector{Tuple{Vector{Float64},Vector{Float64}}}(undef, length(processed_spectra))
-            hw = get(params, :peak_half_window, 10)
-            snr = get(params, :peak_snr, 2.0)
-            intensity_thresh = get(params, :peak_intensity_threshold, 0.0)
-
-            for (i, (mz, y, mode)) in enumerate(processed_spectra)
-                spec_idx = indices[i]
-                if mode == PROFILE
-                    pk_mz, pk_int = detect_peaks_profile(mz, y; half_window=hw, snr_threshold=snr)
-                else # CENTROID
-                    pk_mz, pk_int = detect_peaks_centroid(mz, y; intensity_threshold=intensity_thresh)
-                end
-                peak_results[i] = (pk_mz, pk_int)
-                on_stage(:peaks; idx=spec_idx, mz=pk_mz, intensity=pk_int)
-            end
-            processed_spectra = peak_results # Now contains peak lists
-        
-        elseif step === :align
-            # Alignment requires a reference peak list, typically from the first spectrum
-            reference_peaks = isempty(processed_spectra) ? nothing : processed_spectra[1][1]
-            if reference_peaks === nothing
-                @error "Alignment requires a :peaks step first."; continue
-            end
-            tol = get(params, :align_tolerance, 0.002)
-            for i in 2:length(processed_spectra)
-                tgt_peaks, intens = processed_spectra[i]
-                warp_func = align_peaks_lowess(reference_peaks, tgt_peaks; tolerance=tol)
-                processed_spectra[i] = (warp_func(tgt_peaks), intens)
-                on_stage(:align; idx=indices[i], mz=processed_spectra[i][1], intensity=processed_spectra[i][2])
-            end
-
-        elseif step === :bin
-            all_pks  = [s[1] for s in processed_spectra]
-            all_ints = [s[2] for s in processed_spectra]
-            tol  = get(params, :bin_tolerance, 0.002)
-            freq = get(params, :bin_min_frequency, 0.25)
-            mat, mz_bins = bin_peaks(all_pks, all_ints, tol; frequency_threshold=freq)
-            final_result = FeatureMatrix(mat, mz_bins, indices)
-            break # Binning is the last step
-        end
-    end
-
-    if final_result !== nothing
-        return final_result
-    else
-        @warn "Pipeline finished without a :bin step; returning processed spectra."
-        return processed_spectra
-    end
-end
 
 # =============================================================================
 # 9) Quality Control Metrics
 # =============================================================================
 
 """
-    calculate_ppm_error(measured_mz::Float64, theoretical_mz::Float64) -> Float64
+    get_common_calibration_standards(type::Symbol) -> Dict{Float64, String}
 
-Calculates mass accuracy in parts-per-million (PPM).
+Returns a dictionary of common m/z values for calibration standards.
+
+# Arguments
+- `type::Symbol`: The type of standards to return. Currently supports `:maldi_pos`.
+
+# Returns
+- `Dict{Float64, String}`: A dictionary mapping theoretical m/z to compound names.
+"""
+function get_common_calibration_standards(type::Symbol)
+    if type == :maldi_pos
+        # Common peptide/protein standards for positive mode MALDI
+        return Dict{Float64, String}(
+            1046.5420 => "Angiotensin II",
+            1060.5690 => "Bradykinin",
+            1296.6853 => "Angiotensin I",
+            2465.1989 => "ACTH clip (1-24)"
+        )
+    else
+        @warn "Unsupported calibration standard type: $type. Returning empty dictionary."
+        return Dict{Float64, String}()
+    end
+end
+
+"""
+    calculate_ppm_error(measured_mz::Real, theoretical_mz::Real) -> Float64
+
+Calculates the mass accuracy error in parts-per-million (PPM) between a measured
+and a theoretical m/z value.
+
+# Arguments
+- `measured_mz::Real`: The experimentally measured m/z value.
+- `theoretical_mz::Real`: The known, theoretical m/z value of a compound.
+
+# Returns
+- `Float64`: The calculated PPM error. Returns `Inf` if `theoretical_mz` is zero.
 
 # Formula
-PPM = 10⁶ × |measured_mz - theoretical_mz| / theoretical_mz
+`PPM = 10^6 * |measured_mz - theoretical_mz| / theoretical_mz`
+
+# Example
+```julia
+calculate_ppm_error(100.005, 100.0) # returns 50.0
+```
 """
 function calculate_ppm_error(measured_mz::Real, theoretical_mz::Real)
     if theoretical_mz == 0
@@ -653,30 +651,38 @@ function calculate_ppm_error(measured_mz::Real, theoretical_mz::Real)
 end
 
 """
-    calculate_ppm_error_bulk(measured_mz::Vector{Float64}, theoretical_mz::Vector{Float64}) -> Vector{Float64}
+    calculate_ppm_error_bulk(measured_mz::Vector{<:Real}, theoretical_mz::Vector{<:Real}) -> Vector{Float64}
 
-Calculates PPM errors for multiple mass values.
+Calculates PPM errors for multiple pairs of measured and theoretical mass values.
+
+# Arguments
+- `measured_mz::Vector{<:Real}`: A vector of experimentally measured m/z values.
+- `theoretical_mz::Vector{<:Real}`: A vector of known, theoretical m/z values.
+
+# Returns
+- `Vector{Float64}`: A vector containing the calculated PPM error for each pair.
 """
 function calculate_ppm_error_bulk(measured_mz::Vector{Real}, theoretical_mz::Vector{Real})
     return [calculate_ppm_error(m, t) for (m, t) in zip(measured_mz, theoretical_mz)]
 end
 
 """
-    calculate_resolution_fwhm(mz::Float64, profile_mz::Vector{Float64}, 
-                            profile_intensity::Vector{Float64}) -> Float64
+    calculate_resolution_fwhm(mz::Real, profile_mz::AbstractVector{<:Real}, profile_intensity::AbstractVector{<:Real}) -> Float64
 
-Calculates mass resolution using Full Width at Half Maximum (FWHM).
-
-# Formula
-Resolution = m / Δm, where Δm is FWHM
+Calculates the mass resolution of a peak in profile-mode data using the Full Width at
+Half Maximum (FWHM) method. Resolution is a measure of an instrument's ability to
+distinguish between two peaks of slightly different mass-to-charge ratios.
 
 # Arguments
-- `mz`: Peak centroid m/z
-- `profile_mz`: Full m/z array from profile data
-- `profile_intensity`: Full intensity array from profile data
+- `mz::Real`: The m/z value of the peak's centroid.
+- `profile_mz::AbstractVector{<:Real}`: The full m/z array from the profile-mode spectrum.
+- `profile_intensity::AbstractVector{<:Real}`: The full intensity array from the profile-mode spectrum.
 
 # Returns
-Resolution or NaN if cannot be calculated
+- `Float64`: The calculated resolution (`m / Δm`). Returns `NaN` if the FWHM cannot be determined (e.g., peak is at the edge of the spectrum).
+
+# Formula
+`Resolution = m / Δm`, where `Δm` is the FWHM.
 """
 function calculate_resolution_fwhm(mz::Real, profile_mz::AbstractVector{<:Real},
                                  profile_intensity::AbstractVector{<:Real})
@@ -733,36 +739,48 @@ function find_first_below(v::AbstractVector{<:Real}, threshold::Real)
 end
 
 """
-    analyze_mass_accuracy(msi_data, reference_peaks; ppm_tolerance=20.0)
+    analyze_mass_accuracy(msi_data, reference_peaks; ppm_tolerance=5.0, sample_spectra=100) -> NamedTuple
 
-Analyzes mass accuracy across the dataset using known reference peaks.
+Analyzes the mass accuracy across a sample of spectra from an MSI dataset by comparing
+detected peaks against a list of known reference m/z values. It calculates PPM error
+statistics and suggests an optimal PPM tolerance for future analyses.
 
 # Arguments
-- `msi_data`: Your MSI dataset
-- `reference_peaks`: Dict of theoretical m/z values -> compound names
-- `ppm_tolerance`: Initial tolerance for peak matching
+- `msi_data`: An `MSIData` object containing the dataset.
+- `reference_peaks::Dict{Float64,String}`: A dictionary mapping theoretical m/z values to compound names.
+
+# Keyword Arguments
+- `ppm_tolerance::Float64`: The initial PPM tolerance used to match detected peaks to reference peaks. Defaults to 5.0.
+- `sample_spectra::Int`: The number of spectra to sample from the dataset for the analysis. Defaults to 100.
 
 # Returns
-Comprehensive mass accuracy report
+- `NamedTuple`: A comprehensive report containing statistics like `mean_ppm`, `std_ppm`, `optimal_ppm`, `n_matches`, and lists of matched peaks and all PPM errors.
 """
 function analyze_mass_accuracy(msi_data, reference_peaks::Dict{Float64,String}; 
-                             ppm_tolerance::Float64=5.0, sample_spectra=100)
+                             ppm_tolerance::Float64=5.0, sample_spectra=100,
+                             spectrum_indices=nothing)
     println("\n[ MASS ACCURACY ANALYSIS ]")
     println("PPM tolerance: ", ppm_tolerance, " ppm")
-    println("Spectra to sample: ", sample_spectra)
     
     theoretical_mz = sort(collect(keys(reference_peaks)))
     ppm_errors = Float64[]
     matched_peaks = Tuple{Float64,Float64,String}[]  # (theoretical, measured, compound)
     
-    # Sample spectra across the dataset
-    if sample_spectra >= length(msi_data.spectra_metadata)
-        spectrum_indices = 1:length(msi_data.spectra_metadata)
+    # Determine which spectra to process
+    local indices_to_process
+    if spectrum_indices === nothing
+        println("Spectra to sample: ", sample_spectra)
+        if sample_spectra >= length(msi_data.spectra_metadata)
+            indices_to_process = 1:length(msi_data.spectra_metadata)
+        else
+            indices_to_process = round.(Int, range(1, length(msi_data.spectra_metadata), length=sample_spectra))
+        end
     else
-        spectrum_indices = round.(Int, range(1, length(msi_data.spectra_metadata), length=sample_spectra))
+        indices_to_process = spectrum_indices
+        println("Analyzing $(length(indices_to_process)) specified spectra.")
     end
     
-    for idx in spectrum_indices
+    for idx in indices_to_process
         mz, intensity = GetSpectrum(msi_data, idx)
         
         # Detect peaks in this spectrum
@@ -788,7 +806,7 @@ function analyze_mass_accuracy(msi_data, reference_peaks::Dict{Float64,String};
     
     if isempty(ppm_errors)
         @warn "No peaks matched within $ppm_tolerance ppm tolerance"
-        return (mean_ppm=NaN, std_ppm=NaN, min_ppm=NaN, max_ppm=NaN, optimal_ppm=NaN, n_matches=0, matched_peaks=[], all_ppm_errors=[])
+        return (mean_ppm=NaN, std_ppm=NaN, min_ppm=NaN, max_ppm=NaN, optimal_ppm=NaN, n_matches=0, matched_peaks=[], all_ppm_errors=Float64[])
     end
     
     # Calculate statistics
@@ -813,62 +831,28 @@ function analyze_mass_accuracy(msi_data, reference_peaks::Dict{Float64,String};
 end
 
 """
-    get_common_calibration_standards(standard_type::Symbol)
+    generate_qc_report(msi_data, filename; reference_peaks, output_dir, sample_spectra) -> Tuple
 
-Returns common calibration masses for different instrument types.
+Generates a comprehensive Quality Control (QC) report for an MSI dataset. It assesses
+mass accuracy and resolution, saving the results to CSV files and a summary text file.
 
-# Supported standards
-- `:maldi_pos`: Common MALDI-TOF positive mode calibrants
-- `:maldi_neg`: Common MALDI-TOF negative mode calibrants  
-- `:esi_pos`: ESI positive mode calibrants
-- `:lcms`: LC-MS commonly used standards
+# Arguments
+- `msi_data`: The `MSIData` object to be analyzed.
+- `filename::String`: The name of the input data file, used for reporting.
+
+# Keyword Arguments
+- `reference_peaks::Dict`: A dictionary of reference peaks for mass accuracy analysis. If `nothing`, defaults are loaded using `get_common_calibration_standards`.
+- `output_dir::String`: The directory where the report files will be saved. Defaults to `"qc_results"`.
+- `sample_spectra::Int`: The number of spectra to sample for the analysis. Defaults to 100.
+
+# Returns
+- `Tuple`: A tuple containing the detailed `accuracy_report` (NamedTuple) and `resolution_results` (Vector).
+
+# Side Effects
+- Creates a directory specified by `output_dir`.
+- Writes `mass_accuracy_results.csv`, `resolution_results.csv`, and `qc_summary.txt` into the output directory.
 """
-function get_common_calibration_standards(standard_type::Symbol=:maldi_pos)
-    standards = Dict{Float64,String}()
-    
-    if standard_type == :maldi_pos
-        standards = Dict(
-            104.10754 => "C5H4N2 (Imidazole)",
-            175.11995 => "C6H15O4P (Glycerophosphocholine fragment)",
-            226.15687 => "C10H20NO4P (Phosphocholine)",
-            322.04810 => "[Glu1]-Fibrinopeptide B fragment",
-            379.09247 => "C12H22O11 (Sucrose)",
-            515.32539 => "C26H52NO7P (PC(16:0/0:0))",
-            622.02896 => "C20H12O5S2 (1-Hydroxypyrene-3,6,8-trisulfate)",
-            757.39917 => "C37H74NO8P (PC(34:1))",
-            1046.54198 => "Angiotensin I",
-            1296.68477 => "ACTH clip 1-17",
-            1570.67744 => "ACTH clip 18-39",
-            2465.19829 => "ACTH clip 7-38"
-        )
-    elseif standard_type == :maldi_neg
-        standards = Dict(
-            112.98563 => "C2F3O2 (Trifluoroacetate)",
-            152.99568 => "C2F6S (Perfluoroethylsulfonate)",
-            214.00166 => "C4F7O2 (Heptafluorobutyrate)",
-            264.93278 => "C6F6 (Hexafluorobenzene)",
-            362.96198 => "C8F15O2 (Perfluorooctanoate)",
-            466.96714 => "C10F17O2S (Perfluorooctanesulfonate)"
-        )
-    elseif standard_type == :esi_pos
-        standards = Dict(
-            118.08626 => "C5H12NO2 (Valine)",
-            175.11900 => "C6H15O4P (Phosphocholine fragment)", 
-            524.26496 => "C23H48NO7P (LysoPC(16:0))",
-            622.02896 => "C20H12O5S2 (Standard)",
-            922.00980 => "C18H18O6N3S3 (Ultramark 1621)"
-        )
-    end
-    
-    return standards
-end
-
-"""
-    generate_qc_report(msi_data; reference_peaks, output_dir)
-
-Generates a comprehensive QC report including mass accuracy and resolution.
-"""
-function generate_qc_report(msi_data, filename::String; reference_peaks=nothing, output_dir="qc_results", sample_spectra=100)
+function generate_qc_report(msi_data, filename::String; reference_peaks=nothing, output_dir="qc_results", sample_spectra=100, spectrum_indices=nothing)
     println("\n[ QC REPORT GENERATION ]")
     println("Input file: ", filename)
     println("Output directory: ", output_dir)
@@ -883,7 +867,7 @@ function generate_qc_report(msi_data, filename::String; reference_peaks=nothing,
     println("Using $(length(reference_peaks)) reference masses")
     
     # 1. Analyze mass accuracy
-    accuracy_report = analyze_mass_accuracy(msi_data, reference_peaks, sample_spectra=sample_spectra)
+    accuracy_report = analyze_mass_accuracy(msi_data, reference_peaks; sample_spectra=sample_spectra, spectrum_indices=spectrum_indices)
     
     println("\n" * "="^50)
     println("MASS ACCURACY REPORT")
@@ -976,4 +960,3 @@ function generate_qc_report(msi_data, filename::String; reference_peaks=nothing,
     println("\nQC report saved to: $output_dir")
     return accuracy_report, resolution_results
 end
-
diff --git a/src/PreprocessingVersioning.jl b/src/PreprocessingVersioning.jl
new file mode 100644
index 0000000..aa26537
--- /dev/null
+++ b/src/PreprocessingVersioning.jl
@@ -0,0 +1,446 @@
+# src/PreprocessingVersioning.jl
+
+"""
+This module provides a versioning helper functions in our data preprocessing, module
+inspired by the functionality of the R package MALDIquant. It includes
+functions for quality control, intensity transformation, smoothing, baseline correction,
+normalization, peak picking, alignment, and feature matrix generation.
+
+This module is unfinished and may not end up in the release version.
+"""
+
+# =============================================================================
+# Dependencies
+# =============================================================================
+
+using UUIDs, JLD2, CodecBase # For preprocessing versioning
+
+# =============================================================================
+# Data Structures
+# =============================================================================
+
+"""
+    VersionedSpectralData
+
+A struct to hold a snapshot of spectral data at a specific point in the
+preprocessing workflow. Each transformation creates a new instance of this
+struct, forming a history of all operations.
+
+# Fields
+- `version_id::String`: A unique identifier for this version of the data.
+- `parent_id::String`: The identifier of the version from which this one was derived.
+- `spectra::Vector`: The spectral data itself, typically a vector of `(mz, intensity)` tuples.
+- `processing_step::String`: A description of the operation that created this version (e.g., "smooth", "baseline").
+- `timestamp::DateTime`: The time at which this version was created.
+- `parameters::Dict`: A dictionary of the parameters used in the processing step.
+"""
+struct VersionedSpectralData
+    version_id::String
+    parent_id::String
+    spectra::Vector
+    processing_step::String
+    timestamp::DateTime
+    parameters::Dict
+end
+
+"""
+    MSISession
+
+Manages the state of a preprocessing session, including the history of all
+data versions and the current working version.
+
+# Fields
+- `session_id::String`: A unique identifier for the entire session.
+- `original_data::VersionedSpectralData`: The initial, unprocessed spectral data.
+- `history::Vector{VersionedSpectralData}`: A chronological list of all data versions created during the session.
+- `current_version::VersionedSpectralData`: The version of the data currently being worked on or displayed.
+- `processing_steps::Vector{String}`: A human-readable log of the processing steps applied.
+"""
+mutable struct MSISession
+    session_id::String
+    original_data::VersionedSpectralData
+    history::Vector{VersionedSpectralData}
+    current_version::VersionedSpectralData
+    processing_steps::Vector{String}
+end
+
+# =============================================================================
+# Session Management
+# =============================================================================
+
+"""
+    initialize_session(raw_spectra::Vector, session_name::String) -> MSISession
+
+Creates a new preprocessing session from raw spectral data.
+
+# Arguments
+- `raw_spectra::Vector`: A vector of `(mz, intensity)` tuples representing the initial dataset.
+- `session_name::String`: A name for the session.
+
+# Returns
+- `MSISession`: A new session object initialized with the provided data.
+"""
+function initialize_session(raw_spectra::Vector, session_name::String)::MSISession
+    original = VersionedSpectralData(
+        string(uuid4()),
+        "root",
+        raw_spectra,
+        "Original Data",
+        now(),
+        Dict()
+    )
+    return MSISession(session_name, original, [original], original, [])
+end
+
+"""
+    get_processing_history(session::MSISession) -> Vector{String}
+
+Returns a list of the processing steps applied during the session.
+"""
+get_processing_history(session::MSISession) = [v.processing_step for v in session.history]
+
+"""
+    undo_last_step!(session::MSISession)
+
+Reverts the session to the state before the last processing step was applied.
+This modifies the session in place.
+"""
+function undo_last_step!(session::MSISession)
+    if length(session.history) > 1
+        pop!(session.history)
+        pop!(session.processing_steps)
+        session.current_version = last(session.history)
+    end
+    return session
+end
+
+"""
+    revert_to_version!(session::MSISession, version_id::String)
+
+Reverts the session to a specific version in its history, discarding all
+subsequent changes. This modifies the session in place.
+"""
+function revert_to_version!(session::MSISession, version_id::String)
+    target_idx = findfirst(v -> v.version_id == version_id, session.history)
+    if !isnothing(target_idx)
+        session.history = session.history[1:target_idx]
+        session.current_version = session.history[end]
+        # Also truncate the descriptive processing steps
+        num_steps_to_keep = max(0, target_idx - 1)
+        session.processing_steps = session.processing_steps[1:num_steps_to_keep]
+    end
+    return session
+end
+
+# =============================================================================
+# Versioned Preprocessing
+# =============================================================================
+
+"""
+    apply_processing_step(spectra::Vector, step::Symbol, params::Dict) -> Vector
+
+A dispatcher that applies a single, specified preprocessing function to a set of spectra.
+This is the core function called by the versioned workflow.
+
+# Arguments
+- `spectra::Vector`: The input spectral data.
+- `step::Symbol`: The symbol representing the processing step (e.g., `:smooth`, `:baseline`).
+- `params::Dict`: A dictionary of parameters for the step.
+
+# Returns
+- `Vector`: The processed spectral data.
+"""
+function apply_processing_step(spectra::Vector, step::Symbol, params::Dict)
+    processed = deepcopy(spectra)
+    params_sym = Dict(Symbol(k) => v for (k,v) in params) # Ensure keys are symbols
+
+    if step === :transform
+        for i in eachindex(processed)
+            mz, y = processed[i]
+            processed[i] = (mz, transform_intensity(y; params_sym...))
+        end
+    elseif step === :smooth
+        for i in eachindex(processed)
+            mz, y = processed[i]
+            processed[i] = (mz, smooth_spectrum(y; params_sym...))
+        end
+    elseif step === :baseline
+        for i in eachindex(processed)
+            mz, y = processed[i]
+            baseline = snip_baseline(y; params_sym...)
+            processed[i] = (mz, max.(0.0, y .- baseline))
+        end
+    elseif step === :normalize
+        mode = get(params, :normalize_method, :tic)
+        if mode === :pqn
+            matrix = hcat([float.(p[2]) for p in processed]...)
+            matrix_norm = pqn_normalize(matrix)
+            for i in eachindex(processed)
+                mz, _ = processed[i]
+                processed[i] = (mz, view(matrix_norm, :, i))
+            end
+        else # :tic or :median
+            norm_func = (mode === :tic) ? tic_normalize : median_normalize
+            for i in eachindex(processed)
+                mz, y = processed[i]
+                processed[i] = (mz, norm_func(y))
+            end
+        end
+    elseif step === :peaks
+        peak_results = Vector{Tuple{Vector{Float64},Vector{Float64}}}(undef, length(processed))
+        for (i, (mz, y)) in enumerate(processed)
+            pk_mz, pk_int = detect_peaks_profile(mz, y; params_sym...)
+            peak_results[i] = (pk_mz, pk_int)
+        end
+        return peak_results
+    elseif step === :peaks_wavelet
+        peak_results = Vector{Tuple{Vector{Float64},Vector{Float64}}}(undef, length(processed))
+        for (i, (mz, y)) in enumerate(processed)
+            pk_mz, pk_int = detect_peaks_wavelet(mz, y; params_sym...)
+            peak_results[i] = (pk_mz, pk_int)
+        end
+        return peak_results
+    elseif step === :calibrate
+        return calibrate_spectra(processed; params_sym...)
+    else
+        @warn "Unsupported processing step: $step"
+    end
+    return processed
+end
+
+"""
+    run_versioned_preprocessing(session::MSISession, step::Symbol, params::Dict) -> VersionedSpectralData
+
+Creates a new version of spectral data by applying a single processing step to the
+current version in the session.
+
+# Arguments
+- `session::MSISession`: The current processing session.
+- `step::Symbol`: The processing step to apply.
+- `params::Dict`: Parameters for the processing step.
+
+# Returns
+- `VersionedSpectralData`: A new data version with the transformation applied.
+"""
+function run_versioned_preprocessing(session::MSISession, step::Symbol, params::Dict)
+    # Create a new version based on the current one
+    parent_version = session.current_version
+    new_version_id = string(uuid4())
+    
+    # Apply the processing step
+    new_spectra = apply_processing_step(parent_version.spectra, step, params)
+    
+    # Create the new versioned data object
+    new_version = VersionedSpectralData(
+        new_version_id,
+        parent_version.version_id,
+        new_spectra,
+        string(step),
+        now(),
+        params
+    )
+    
+    return new_version
+end
+
+"""
+    apply_processing_with_version!(session::MSISession, step::Symbol, ui_params::Dict)
+
+A high-level function to apply a processing step, create a new version, and
+update the session state in place.
+
+# Arguments
+- `session::MSISession`: The session to modify.
+- `step::Symbol`: The processing step to apply.
+- `ui_params::Dict`: A dictionary of parameters, typically from a UI.
+"""
+function apply_processing_with_version!(session::MSISession, step::Symbol, ui_params::Dict)
+    # In a real app, you would validate and convert UI params here.
+    processing_params = ui_params
+    
+    # Create the new version
+    new_version = run_versioned_preprocessing(session, step, processing_params)
+    
+    # Update the session
+    push!(session.history, new_version)
+    session.current_version = new_version
+    
+    # Describe the step for the history log
+    step_description = "$(string(step)) with params: " * join(["$k=$v" for (k,v) in processing_params], ", ")
+    push!(session.processing_steps, step_description)
+    
+    return session
+end
+
+# =============================================================================
+# Data Serialization & Export
+# =============================================================================
+
+"""
+    save_spectral_version(data::VersionedSpectralData, filepath::String)
+
+Saves a `VersionedSpectralData` object to a file using the JLD2 format.
+"""
+function save_spectral_version(data::VersionedSpectralData, filepath::String)
+    JLD2.save_object(filepath, data)
+end
+
+"""
+    load_spectral_version(filepath::String) -> VersionedSpectralData
+
+Loads a `VersionedSpectralData` object from a JLD2 file.
+"""
+function load_spectral_version(filepath::String)::VersionedSpectralData
+    return JLD2.load_object(filepath)
+end
+
+"""
+    _encode_base64(data::AbstractVector{<:Real}) -> String
+
+Helper function to convert a numeric vector into a Base64 encoded string.
+"""
+function _encode_base64(data::AbstractVector{T}) where T <: Real
+    bytes = reinterpret(UInt8, data)
+    return base64encode(bytes)
+end
+
+"""
+    export_to_mzml(session::MSISession, filepath::String)
+
+Exports the current version of the spectral data in a session to a standard
+.mzML file.
+
+# Arguments
+- `session::MSISession`: The current session.
+- `filepath::String`: The path for the output .mzML file.
+"""
+function export_to_mzml(session::MSISession, filepath::String)
+    spectra_to_export = session.current_version.spectra
+    
+    open(filepath, "w") do f
+        # XML Header
+        write(f, "<?xml version=\"1.0\" encoding=\"UTF-8\"?>\n")
+        write(f, "<mzML xmlns=\"http://psi.hupo.org/ms/mzml\" xmlns:xsi=\"http://www.w3.org/2001/XMLSchema-instance\" xsi:schemaLocation=\"http://psi.hupo.org/ms/mzml http://psidev.info/files/ms/mzML/xsd/mzML1.1.0.xsd\" version=\"1.1\">\n")
+        
+        # CV List
+        write(f, "  <cvList count=\"2\">\n")
+        write(f, "    <cv id=\"MS\" fullName=\"Proteomics Standards Initiative Mass Spectrometry Ontology\" version=\"4.1.123\" URI=\"https://raw.githubusercontent.com/HUPO-PSI/psi-ms-CV/master/psi-ms.obo\"/>\n")
+        write(f, "    <cv id=\"UO\" fullName=\"Unit Ontology\" version=\"07:03:2022\" URI=\"https://raw.githubusercontent.com/bio-ontology-research-group/unit-ontology/master/unit.obo\"/>\n")
+        write(f, "  </cvList>\n")
+
+        # File Description
+        write(f, "  <fileDescription>\n    <fileContent>\n")
+        write(f, "      <cvParam cvRef=\"MS\" accession=\"MS:1000579\" name=\"MS1 spectrum\"/>\n")
+        write(f, "    </fileContent>\n  </fileDescription>\n")
+
+        # Run and Spectrum List
+        write(f, "  <run id=\"run1\" defaultInstrumentConfigurationRef=\"instrument1\">\n")
+        write(f, "    <spectrumList count=\"$(length(spectra_to_export))\" defaultDataProcessingRef=\"dp1\">\n")
+
+        for (i, (mz, intensity)) in enumerate(spectra_to_export)
+            # Ensure data is in the correct format for encoding
+            mz_64 = convert(Vector{Float64}, mz)
+            int_32 = convert(Vector{Float32}, intensity)
+
+            # Base64 encode the binary data
+            mz_b64 = _encode_base64(mz_64)
+            int_b64 = _encode_base64(int_32)
+
+            write(f, "      <spectrum index=\"$(i-1)\" id=\"scan=$(i)\" defaultArrayLength=\"$(length(mz))\">\n")
+            write(f, "        <cvParam cvRef=\"MS\" accession=\"MS:1000511\" name=\"ms level\" value=\"1\"/>\n")
+            # Assuming profile mode for processed data, could be made dynamic
+            write(f, "        <cvParam cvRef=\"MS\" accession=\"MS:1000128\" name=\"profile spectrum\"/>\n")
+            
+            write(f, "        <binaryDataArrayList count=\"2\">\n")
+            
+            # m/z array
+            write(f, "          <binaryDataArray encodedLength=\"$(length(mz_b64))\">\n")
+            write(f, "            <cvParam cvRef=\"MS\" accession=\"MS:1000514\" name=\"m/z array\" unitCvRef=\"MS\" unitAccession=\"MS:1000040\" unitName=\"m/z\"/>\n")
+            write(f, "            <cvParam cvRef=\"MS\" accession=\"MS:1000523\" name=\"64-bit float\"/>\n")
+            write(f, "            <cvParam cvRef=\"MS\" accession=\"MS:1000576\" name=\"no compression\"/>\n")
+            write(f, "            <binary>$(mz_b64)</binary>\n")
+            write(f, "          </binaryDataArray>\n")
+
+            # Intensity array
+            write(f, "          <binaryDataArray encodedLength=\"$(length(int_b64))\">\n")
+            write(f, "            <cvParam cvRef=\"MS\" accession=\"MS:1000515\" name=\"intensity array\" unitCvRef=\"MS\" unitAccession=\"MS:1000131\" unitName=\"number of detector counts\"/>\n")
+            write(f, "            <cvParam cvRef=\"MS\" accession=\"MS:1000521\" name=\"32-bit float\"/>\n")
+            write(f, "            <cvParam cvRef=\"MS\" accession=\"MS:1000576\" name=\"no compression\"/>\n")
+            write(f, "            <binary>$(int_b64)</binary>\n")
+            write(f, "          </binaryDataArray>\n")
+            
+            write(f, "        </binaryDataArrayList>\n")
+            write(f, "      </spectrum>\n")
+        end
+
+        write(f, "    </spectrumList>\n")
+        write(f, "  </run>\n")
+        
+        # Dummy instrument and data processing info
+        write(f, "  <instrumentConfigurationList count=\"1\">\n")
+        write(f, "    <instrumentConfiguration id=\"instrument1\"/>\n")
+        write(f, "  </instrumentConfigurationList>\n")
+        write(f, "  <dataProcessingList count=\"1\">\n")
+        write(f, "    <dataProcessing id=\"dp1\"/>\n")
+        write(f, "  </dataProcessingList>\n")
+
+        write(f, "</mzML>\n")
+    end
+    println("Successfully exported current data to $filepath")
+end
+
+
+# =============================================================================
+# Stubs for UI Integration
+# =============================================================================
+
+"""
+    validate_ui_parameters(step::Symbol, ui_params::Dict) -> Tuple{Bool, String}
+
+(Stub) Validates parameters from a UI before they are used in a processing step.
+"""
+function validate_ui_parameters(step::Symbol, ui_params::Dict)::Tuple{Bool, String}
+    # In a real implementation, this would check types, ranges, etc.
+    # For example, for :smooth, ensure 'sg_window' is an odd integer.
+    println("Validating parameters for step: $step")
+    return (true, "Parameters are valid.")
+end
+
+"""
+    run_processing_with_progress(session, step, params, progress_callback)
+
+(Stub) A wrapper for running a processing step that includes a progress reporting callback.
+"""
+function run_processing_with_progress(session, step, params, progress_callback)
+    progress_callback(0.0, "Starting $step...")
+    
+    # This is a simplified example. Real implementation would need to
+    # hook into the loops inside apply_processing_step.
+    new_version = run_versioned_preprocessing(session, step, params)
+    
+    progress_callback(1.0, "Finished $step.")
+    return new_version
+end
+
+"""
+    safe_processing_application!(session::MSISession, step::Symbol, params::Dict)
+
+(Stub) A safe wrapper to apply a processing step that includes error handling and rollback.
+"""
+function safe_processing_application!(session::MSISession, step::Symbol, params::Dict)
+    num_history = length(session.history)
+    try
+        println("Safely applying step: $step")
+        apply_processing_with_version!(session, step, params)
+    catch e
+        @error "Processing step '$step' failed!" exception=(e, catch_backtrace())
+        # Roll back to the previous state
+        while length(session.history) > num_history
+            undo_last_step!(session)
+        end
+        println("Session has been rolled back to the previous state.")
+        return session # Return the rolled-back session
+    end
+    return session # Return the updated session
+end
+
diff --git a/test/readme.md b/test/readme.md
index 197af0d..1965e10 100644
--- a/test/readme.md
+++ b/test/readme.md
@@ -40,10 +40,10 @@ Before running the tests, you must edit the `test/run_tests.jl` file to point to
 2.  **Execute the Different Test Scripts**:
     Run the following command from the project's root directory. This will install the necessary dependencies and run the tests.
     ```bash
-    julia --project=. test/run_tests.jl
+    julia --threads auto --project=. test/run_tests.jl
     ```
     ```bash
-    julia --project=. test/run_preprocessing.jl
+    julia --threads auto --project=. test/run_preprocessing.jl
     ```
 
 3.  **Check the Results**:
diff --git a/test/run_preprocessing.jl b/test/run_preprocessing.jl
index 4ea526a..bc70023 100644
--- a/test/run_preprocessing.jl
+++ b/test/run_preprocessing.jl
@@ -1,468 +1,371 @@
 # test/run_preprocessing.jl
 
 # ===================================================================
-# Test Environment for the Preprocessing.jl Module
+# Preprocessing Test Suite for JuliaMSI
 # ===================================================================
-# This script tests the full preprocessing pipeline on single spectra
-# and total spectra from both .mzML and .imzML files.
-# It generates an overlay plot showing all preprocessing stages and
-# saves the resulting feature matrix to a CSV file.
+# This script provides a customizable workflow to test and visualize
+# the effects of different preprocessing steps on individual spectra
+# from .mzML or .imzML files.
 #
 # Instructions:
-# 1. Ensure the file paths in the "CONFIG" section are correct.
+# 1. Configure the file paths and parameters in the "CONFIG" section.
 # 2. Run the script from the project's root directory:
-#    julia test/run_preprocessing.jl
-# 3. Check the `test/results/` folder for output plots and CSVs.
+#    julia --threads auto --project=. test/run_preprocessing.jl
+# 3. Check the configured `RESULTS_DIR` for output plots and CSV files.
 # ===================================================================
 
 using Printf
 using CairoMakie
+using DataFrames
+using CSV
+using Statistics
+using Interpolations
 import Pkg
-using DataFrames # For saving FeatureMatrix to CSV
-using CSV        # For saving FeatureMatrix to CSV
-using Statistics # For mean()
 
-# --- Load Modules ---
-# Activate the project environment to access dependencies
+# --- Load the MSI_src Module ---
+# Activate the project environment at the parent directory of this test script
 Pkg.activate(joinpath(@__DIR__, ".."))
-using MSI_src # This brings in Preprocessing.jl functions via export
+using MSI_src
+
 
 # ===================================================================
-# CONFIG: Test files and parameters
+# CONFIG: CUSTOMIZE YOUR PREPROCESSING WORKFLOW HERE
 # ===================================================================
 
-# --- Test Files ---
-# An mzML file for testing spectrum-based processing
-# const TEST_MZML_FILE = "/home/pixel/Documents/Cinvestav_2025/Analisis/CE4_BF_R1/CE4_BF_R1.mzML"
-# const TEST_MZML_FILE = "/home/pixel/Documents/Cinvestav_2025/Analisis/set de datos MS/Leaf_profile_LD_LTP_MS.mzML"
-const TEST_MZML_FILE = "/home/pixel/Documents/Cinvestav_2025/Analisis/set de datos MS/Escopolamina_tuneo_fraq_20ev.mzML"
-#const TEST_MZML_FILE = "/home/pixel/Documents/Cinvestav_2025/Analisis/set de datos MS/Atropina_tuneo_fraq_20ev.mzML"
+# --- Input and Output ---
+const TEST_FILE = "/home/pixel/Documents/Cinvestav_2025/Analisis/salida/Stomach_DHB_uncompressed.imzML" # Can be .mzML or .imzML
+const RESULTS_DIR = "test/results/preprocessing"
+const NUM_SPECTRA_TO_PROCESS = nothing # Set to `nothing` to process all spectra
 
-const MZML_SPECTRUM_ID = 1
+# --- Internal Standards for Calibration and QC ---
+# replace these with m/z values of known compounds present in your dataset.
+const INTERNAL_STANDARDS = Dict(
+    "P13" => 432.6584,
+    "P15" => 464.6059,
+    "P17" => 526.5534,
+    "P21" => 650.4485,
+    "P25" => 774.3435,
+    "P29" => 898.2385,
+    "P31" => 950.1861,
+    "P33" => 1022.1336,
+    "P37" => 1146.0286,
+    "P45" => 1593.8187,
+    "Unknown 1" => 772.433,
+    "Unknown 2" => 772.5253
+)
 
-# An imzML file for testing
-# const TEST_IMZML_FILE = "/home/pixel/Documents/Cinvestav_2025/Analisis/CE4_BF_R1/CE4_BF_R1.imzML"
-# const IMZML_COORDS = (50, 50) 
-const TEST_IMZML_FILE = "/home/pixel/Documents/Cinvestav_2025/Analisis/salida/Stomach_DHB_uncompressed.imzML"
-const IMZML_COORDS = (1997, 639) 
+# --- Preprocessing Step Parameters ---
 
-# --- Output Directory ---
-const RESULTS_DIR = "test/results"
+# Step 0: Quality Control (QC)
+const QC_PARAMS = (
+    ppm_tolerance = 5.0, # PPM tolerance for matching internal standards
+)
+
+# Step 1: Calibration
+const CALIBRATION_PARAMS = (
+    enabled = true,
+    ppm_tolerance = 5.0, # PPM tolerance for finding calibration peaks
+)
+
+# Step 2: Smoothing
+# Note: Smoothing is less effective and often unnecessary for centroid data.
+const SMOOTHING_PARAMS = (
+    enabled = false,
+    window = 9,
+    order = 2,
+)
+
+# Step 3: Baseline Correction
+# Note: Baseline correction is less effective and often unnecessary for centroid data.
+const BASELINE_CORRECTION_PARAMS = (
+    enabled = false,
+    iterations = 100,
+)
+
+# Step 4: Normalization
+const NORMALIZATION_PARAMS = (
+    enabled = true,
+    method = :tic, # :tic, :median, or :none
+)
+
+# Step 5: Peak Detection
+const PEAK_DETECTION_PARAMS = (
+    enabled = true,
+    method = :centroid, # Use :centroid for centroided data, :profile for profile data
+    # --- Parameters for :centroid method ---
+    intensity_threshold = 0.0, # Filters out peaks below this absolute intensity
+    # --- Parameters for :profile method ---
+    half_window = 10,
+    snr_threshold = 3.0,
+    min_peak_prominence = 0.05,
+    merge_peaks_tolerance = 0.002,
+)
+
+# Step 6: Peak Alignment (Warping)
+# This is performed after collecting peaks from all spectra.
+const PEAK_ALIGNMENT_PARAMS = (
+    enabled = true,
+    tolerance = 0.002,
+    tolerance_unit = :mz, # :mz or :ppm
+    min_matched_peaks = 3, # Lowered for sparse centroid data
+)
+
+# Step 7: Peak Binning
+const PEAK_BINNING_PARAMS = (
+    enabled = true,
+    tolerance = 0.1,
+    tolerance_unit = :mz, # :mz or :ppm
+    frequency_threshold = 0.1, # Min fraction of spectra a peak must be in
+)
 
 # ===================================================================
-# HELPER FUNCTIONS FOR PLOTTING
+# HELPER FUNCTIONS
 # ===================================================================
 
 """
-    plot_overlay_stages(collected_data, output_path, title)
+    plot_spectrum_on_fig(fig_pos, mz, intensity, title)
 
-Creates a single plot overlaying spectra from different preprocessing stages.
+Helper function to plot a spectrum on a specific position of a CairoMakie figure.
 """
-function plot_overlay_stages(collected_data, output_path, title)
-    fig = Figure(size = (1400, 800))
-    ax = Axis(fig[1, 1], title=title, xlabel="m/z", ylabel="Intensity")
-
-    colors = Makie.wong_colors() # A good set of distinct colors
-
-    for (i, (stage, mz, intensity)) in enumerate(collected_data)
-        color = colors[mod1(i, length(colors))] # Cycle through colors
-        
-        # Plot the spectrum as a line
-        lines!(ax, mz, intensity, color=color, label=string(stage))
-
-        # If it's the peaks stage, also mark the peak tops
-        if stage == :peaks
-            scatter!(ax, mz, intensity, color=color, marker=:circle, markersize=8, label="$(string(stage)) (tops)")
-        end
-    end
-    axislegend(ax, position=:rt) # Right top position
-    save(output_path, fig)
-    println("SUCCESS: Overlay plot saved to $output_path")
+function plot_spectrum_on_fig(fig_pos, mz, intensity, title)
+    ax = Axis(fig_pos, title=title, xlabel="m/z", ylabel="Intensity")
+    lines!(ax, mz, intensity)
 end
 
 
 # ===================================================================
-# TEST DEFINITIONS
+# MAIN PROCESSING SCRIPT
 # ===================================================================
 
-"""
-    test_full_pipeline(msi_data, spectrum_id; output_dir, file_type_prefix)
+function run_preprocessing_suite()
+    println("Starting Preprocessing Test Suite...")
+    mkpath(RESULTS_DIR)
 
-Tests the full preprocessing pipeline on a single spectrum and saves a plot
-for each intermediate step using the `on_stage` callback.
-`spectrum_id` can be an `Int` (for mzML) or a `Tuple{Int, Int}` (for imzML).
-"""
-function test_full_pipeline(msi_data, spectrum_id; output_dir, file_type_prefix, mz_tolerance=0.002)
-    println("\n--- Testing Full Preprocessing Pipeline on Spectrum: $spectrum_id (File Type: $file_type_prefix) ---")
+    # --- Load Data ---
+    println("Loading data from: $TEST_FILE")
+    if !isfile(TEST_FILE)
+        println("ERROR: Test file not found. Please check the path in the CONFIG section.")
+        return
+    end
+    msi_data = @time OpenMSIData(TEST_FILE)
+    println("Data loaded successfully. Found $(length(msi_data.spectra_metadata)) spectra.")
 
-    # 1. Determine the spectrum index
-    local spec_idx
-    if spectrum_id isa Int
-        spec_idx = spectrum_id
-    else # Tuple for imzML
-        spec_idx = msi_data.coordinate_map[spectrum_id...]
+    # --- Determine which spectra to process ---
+    total_spectra = length(msi_data.spectra_metadata)
+    indices_to_process = if NUM_SPECTRA_TO_PROCESS === nothing
+        1:total_spectra
+    else
+        unique(round.(Int, range(1, total_spectra, length=min(NUM_SPECTRA_TO_PROCESS, total_spectra))))
+    end
+    println("Will process $(length(indices_to_process)) spectra.")
+
+    # --- Data storage for results ---
+    qc_results = DataFrame(spectrum_idx=Int[], metric=String[], value=Float64[], compound=String[])
+    all_processed_peaks_mz = Vector{Vector{Float64}}()
+    all_processed_peaks_int = Vector{Vector{Float64}}()
+    spectrum_indices_with_peaks = Int[]
+
+    # --- Main Loop: Process each spectrum individually ---
+    println("\n" * "="^20 * " Processing Individual Spectra " * "="^20)
+    for (i, idx) in enumerate(indices_to_process)
+        print("\rProcessing spectrum #$idx ($(i)/$(length(indices_to_process)))...")
+        
+        process_spectrum(msi_data, idx) do mz, intensity
+            if qc_is_empty(mz, intensity) || !qc_is_regular(mz)
+                # @warn "Skipping empty or irregular spectrum #$idx"
+                return
+            end
+
+            original_mz, original_intensity = copy(mz), copy(intensity)
+            processed_mz, processed_intensity = copy(mz), copy(intensity)
+
+            # --- Step 0: QC (PPM and Resolution) ---
+            # This happens before any modification to the m/z axis
+            let ref_masses = collect(values(INTERNAL_STANDARDS))
+                matched_peaks = find_calibration_peaks(processed_mz, processed_intensity, ref_masses, ppm_tolerance=QC_PARAMS.ppm_tolerance)
+                
+                for (compound, theoretical_mz) in INTERNAL_STANDARDS
+                    # Find the measured m/z that corresponds to this theoretical_mz
+                    measured_mz = 0.0
+                    for (theo, meas) in matched_peaks
+                        if theo == theoretical_mz
+                            measured_mz = meas
+                            break
+                        end
+                    end
+
+                    if measured_mz > 0
+                        # PPM Error
+                        ppm_error = calculate_ppm_error(measured_mz, theoretical_mz)
+                        push!(qc_results, (idx, "ppm_error", ppm_error, compound))
+                        
+                        # Resolution
+                        resolution = calculate_resolution_fwhm(measured_mz, processed_mz, processed_intensity)
+                        if !isnan(resolution)
+                            push!(qc_results, (idx, "resolution", resolution, compound))
+                        end
+                    end
+                end
+            end
+
+            # --- Step 1: Calibration ---
+            if CALIBRATION_PARAMS.enabled
+                ref_masses = collect(values(INTERNAL_STANDARDS))
+                matched_peaks = find_calibration_peaks(processed_mz, processed_intensity, ref_masses, ppm_tolerance=CALIBRATION_PARAMS.ppm_tolerance)
+                if length(matched_peaks) >= 2
+                    measured = sort(collect(values(matched_peaks)))
+                    theoretical = sort(collect(keys(matched_peaks)))
+                    itp = linear_interpolation(measured, theoretical, extrapolation_bc=Line())
+                    processed_mz = itp(processed_mz)
+                end
+            end
+
+            # --- Step 2: Smoothing ---
+            if SMOOTHING_PARAMS.enabled
+                processed_intensity = smooth_spectrum(processed_intensity, window=SMOOTHING_PARAMS.window, order=SMOOTHING_PARAMS.order)
+            end
+
+            # --- Step 3: Baseline Correction ---
+            if BASELINE_CORRECTION_PARAMS.enabled
+                baseline = snip_baseline(processed_intensity, iterations=BASELINE_CORRECTION_PARAMS.iterations)
+                processed_intensity .-= baseline
+                processed_intensity = max.(0, processed_intensity) # Ensure non-negativity
+            end
+
+            # --- Step 4: Normalization ---
+            if NORMALIZATION_PARAMS.enabled && NORMALIZATION_PARAMS.method != :none
+                if NORMALIZATION_PARAMS.method == :tic
+                    processed_intensity = tic_normalize(processed_intensity)
+                elseif NORMALIZATION_PARAMS.method == :median
+                    processed_intensity = median_normalize(processed_intensity)
+                end
+            end
+
+            # --- Step 5: Peak Detection ---
+            local pk_mz, pk_int
+            pk_mz, pk_int = Float64[], Float64[] # Initialize empty
+            if PEAK_DETECTION_PARAMS.enabled
+                if PEAK_DETECTION_PARAMS.method == :profile
+                    pk_mz, pk_int = detect_peaks_profile(processed_mz, processed_intensity,
+                        half_window=PEAK_DETECTION_PARAMS.half_window,
+                        snr_threshold=PEAK_DETECTION_PARAMS.snr_threshold,
+                        min_peak_prominence=PEAK_DETECTION_PARAMS.min_peak_prominence,
+                        merge_peaks_tolerance=PEAK_DETECTION_PARAMS.merge_peaks_tolerance)
+                elseif PEAK_DETECTION_PARAMS.method == :wavelet
+                    pk_mz, pk_int = detect_peaks_wavelet(processed_mz, processed_intensity)
+                else # :centroid
+                    pk_mz, pk_int = detect_peaks_centroid(processed_mz, processed_intensity)
+                end
+                
+                if !isempty(pk_mz)
+                    push!(all_processed_peaks_mz, pk_mz)
+                    push!(all_processed_peaks_int, pk_int)
+                    push!(spectrum_indices_with_peaks, idx)
+                end
+            end
+
+            # --- Visualization of a sample spectrum ---
+            if i == 1 # Only plot the first processed spectrum
+                println("\nGenerating example plots for spectrum #$idx...")
+                fig = Figure(size=(1200, 800))
+                
+                plot_spectrum_on_fig(fig[1,1], original_mz, original_intensity, "1. Original Spectrum")
+                plot_spectrum_on_fig(fig[2,1], processed_mz, processed_intensity, "2. After All Steps (Before Peak Picking)")
+                
+                # Plot detected peaks
+                ax = Axis(fig[3,1], title="3. Detected Peaks")
+                if !isempty(pk_mz)
+                    stem!(ax, pk_mz, pk_int)
+                end
+                
+                save(joinpath(RESULTS_DIR, "example_spectrum_processing.png"), fig)
+                println("Saved example processing plots.")
+            end
+        end
+    end
+    println("\nIndividual spectrum processing complete.")
+
+    # --- Save QC Results ---
+    if !isempty(qc_results)
+        println("\n" * "="^20 * " Generating QC Report " * "="^20)
+        CSV.write(joinpath(RESULTS_DIR, "qc_results.csv"), qc_results)
+        println("QC results saved to qc_results.csv")
+
+        # Generate summary plots for QC
+        fig = Figure(size=(1200, 600))
+        
+        # PPM Error Histogram
+        ppm_errors = filter(row -> row.metric == "ppm_error", qc_results).value
+        if !isempty(ppm_errors)
+            ax1 = Axis(fig[1,1], title="PPM Error Distribution", xlabel="PPM Error")
+            hist!(ax1, ppm_errors, bins=30)
+        end
+
+        # Resolution Histogram
+        resolutions = filter(row -> row.metric == "resolution", qc_results).value
+        if !isempty(resolutions)
+            ax2 = Axis(fig[1,2], title="Resolution Distribution", xlabel="Resolution (FWHM)")
+            hist!(ax2, resolutions, bins=30)
+        end
+        
+        save(joinpath(RESULTS_DIR, "qc_summary_plots.png"), fig)
+        println("QC summary plots saved.")
     end
 
-    if spec_idx == 0
-        println("SKIPPED: No spectrum found at coordinates $spectrum_id.")
+    if isempty(all_processed_peaks_mz)
+        @warn "No peaks were detected in any of the processed spectra. Skipping alignment and binning."
         return
     end
 
-    # 2. Define the pipeline steps in the desired order
-    pipeline_steps = [
-        :qc,
-        :transform,
-        :smooth,
-        :baseline,
-        :normalize,
-        :peaks,
-        :align, # Align requires multiple spectra, but we'll run it on a single one for now (will warn)
-        :bin
-    ]
+    # --- Step 6: Peak Alignment ---
+    println("\n" * "="^20 * " Aligning Peaks " * "="^20)
+    aligned_peaks_mz = copy(all_processed_peaks_mz)
+    if PEAK_ALIGNMENT_PARAMS.enabled
+        # Create a reference peak list (e.g., from the spectrum with the most peaks)
+        ref_idx = argmax(length.(all_processed_peaks_mz))
+        reference_peaks = all_processed_peaks_mz[ref_idx]
+        println("Using spectrum $(spectrum_indices_with_peaks[ref_idx]) as alignment reference.")
 
-    # Define parameters for each step
-    params = Dict(
-        :transform_method => :sqrt,
-        :sg_window => 15,
-        :sg_order => 2,
-        :snip_iterations => 100,
-        :normalize_method => :tic,
-        :peak_half_window => 10,
-        :peak_snr => 3.0,
-        :peak_intensity_threshold => 0.0, # For centroid peak detection
-        :align_tolerance => mz_tolerance,
-        :bin_tolerance => mz_tolerance,
-        :bin_min_frequency => 0.0 # Keep all bins for a single spectrum
-    )
+        for i in 1:length(aligned_peaks_mz)
+            if i == ref_idx continue end
+            alignment_func = align_peaks_lowess(reference_peaks, aligned_peaks_mz[i],
+                tolerance=PEAK_ALIGNMENT_PARAMS.tolerance,
+                tolerance_unit=PEAK_ALIGNMENT_PARAMS.tolerance_unit,
+                min_matched_peaks=PEAK_ALIGNMENT_PARAMS.min_matched_peaks)
+            
+            aligned_peaks_mz[i] = alignment_func(aligned_peaks_mz[i])
+        end
+        println("Peak alignment complete.")
+    end
 
-    # 3. Define the on_stage callback to collect data for overlay plot and save separate plots
-    collected_stage_data = []
-    stage_counter = Ref(0) # Initialize counter for sequential naming
-    normalized_spectrum = nothing # Variable to hold the normalized spectrum
-
-    function stage_callback(stage; idx, mz, intensity)
-        stage_counter[] += 1 # Increment counter
-        println("  -> Generating plot for stage: $stage")
+    # --- Step 7: Peak Binning & Feature Matrix Generation ---
+    println("\n" * "="^20 * " Binning Peaks " * "="^20)
+    if PEAK_BINNING_PARAMS.enabled
+        feature_matrix, mz_bins = bin_peaks(aligned_peaks_mz, all_processed_peaks_int,
+            PEAK_BINNING_PARAMS.tolerance,
+            tolerance_unit=PEAK_BINNING_PARAMS.tolerance_unit,
+            frequency_threshold=PEAK_BINNING_PARAMS.frequency_threshold)
         
-        local fig # Make fig available in the whole function scope
-
-        if stage == :normalize
-            normalized_spectrum = (mz, intensity)
-            fig = plot_stage_spectrum(mz, intensity, title="Stage: $stage (Spectrum $spectrum_id)")
-        elseif stage == :peaks && normalized_spectrum !== nothing
-            # For the peaks stage, plot the normalized spectrum as a base layer
-            fig = Figure(size = (1400, 500))
-            ax = Axis(fig[1, 1], title="Stage: Peaks (Spectrum $spectrum_id)", xlabel="m/z", ylabel="Intensity")
-            lines!(ax, normalized_spectrum[1], normalized_spectrum[2], color=:gray, label="Normalized Spectrum")
-            scatter!(ax, mz, intensity, color=:red, marker=:circle, markersize=8, label="Detected Peaks")
-            axislegend(ax)
+        if !isempty(feature_matrix)
+            df = DataFrame(feature_matrix, :auto)
+            rename!(df, ["bin_$(i)" for i in 1:size(df, 2)])
+            insertcols!(df, 1, :spectrum_idx => spectrum_indices_with_peaks)
+            
+            CSV.write(joinpath(RESULTS_DIR, "feature_matrix.csv"), df)
+            println("Feature matrix saved to feature_matrix.csv")
+            
+            # Save bin m/z ranges
+            bin_df = DataFrame(bin_index=1:length(mz_bins), mz_start=[b[1] for b in mz_bins], mz_end=[b[2] for b in mz_bins])
+            CSV.write(joinpath(RESULTS_DIR, "bin_definitions.csv"), bin_df)
+            println("Bin definitions saved to bin_definitions.csv")
         else
-            # Default plotting for all other stages
-            fig = plot_stage_spectrum(mz, intensity, title="Stage: $stage (Spectrum $spectrum_id)")
+            @warn "Feature matrix was empty after binning."
         end
-
-        # Save the figure
-        stage_output_path = joinpath(output_dir, "$(file_type_prefix)_$(spectrum_id)_$(stage_counter[])_$(stage).png")
-        save(stage_output_path, fig)
-
-        # Collect data for overlay plot
-        push!(collected_stage_data, (stage, mz, intensity))
     end
 
-    # 4. Run the pipeline on the single spectrum
-    println("Running pipeline with steps: $pipeline_steps")
-    processed_result = run_preprocessing_pipeline(
-        msi_data, 
-        [spec_idx], # The pipeline expects a vector of indices
-        steps=pipeline_steps,
-        params=params,
-        on_stage=stage_callback
-    )
-
-    # 5. Generate and save the overlay plot
-    overlay_output_path = joinpath(output_dir, "$(file_type_prefix)_$(spectrum_id)_all_stages_overlay.png")
-    plot_overlay_stages(collected_stage_data, overlay_output_path, "Preprocessing Stages Overlay (Spectrum $spectrum_id)")
-
-    # 6. Save feature matrix if generated
-    if processed_result isa FeatureMatrix
-        feature_matrix_output_path = joinpath(output_dir, "$(file_type_prefix)_$(spectrum_id)_feature_matrix.csv")
-        # Convert mz_bins to a more readable format for CSV
-        mz_labels = ["$(round(b[1], digits=4))_$(round(b[2], digits=4))" for b in processed_result.mz_bins]
-        df = DataFrame(processed_result.matrix, Symbol.(mz_labels))
-        CSV.write(feature_matrix_output_path, df)
-        println("SUCCESS: Feature matrix saved to $feature_matrix_output_path")
-    else
-        @warn "Pipeline did not return a FeatureMatrix for Spectrum $spectrum_id."
-        processed_result
-    end
-
-    println("--- Pipeline test finished for Spectrum: $spectrum_id (File Type: $file_type_prefix) ---")
-    println("Check the '$(output_dir)' directory for output plots and CSVs.")
+    println("\nPreprocessing Test Suite finished successfully!")
 end
 
-"""
-    test_full_pipeline_on_total_spectrum(msi_data; output_dir, file_type_prefix)
-
-Tests the full preprocessing pipeline on the *total spectrum* (sum of all spectra)
-and saves a plot for each intermediate step.
-"""
-function test_full_pipeline_on_total_spectrum(msi_data; output_dir, file_type_prefix, mz_tolerance=0.002)
-    println("\n--- Testing Full Preprocessing Pipeline on TOTAL Spectrum (File Type: $file_type_prefix) ---")
-
-    # 1. Get the total spectrum
-    total_mz, total_intensity = get_total_spectrum(msi_data)
-    total_spectrum = (total_mz, total_intensity)
-
-    if qc_is_empty(total_mz, total_intensity)
-        println("SKIPPED: Total spectrum is empty.")
-        return
-    end
-
-    # 2. Define the pipeline steps and parameters (same as for single spectrum)
-    pipeline_steps = [
-        :qc,
-        :transform,
-        :smooth,
-        :baseline,
-        :normalize,
-        :peaks,
-        :align, # Align requires multiple spectra, but we'll run it on a single one for now (will warn)
-        :bin
-    ]
-
-    params = Dict(
-        :transform_method => :sqrt,
-        :sg_window => 15,
-        :sg_order => 2,
-        :snip_iterations => 100,
-        :normalize_method => :tic,
-        :peak_half_window => 10,
-        :peak_snr => 3.0,
-        :align_tolerance => mz_tolerance,
-        :bin_tolerance => mz_tolerance,
-        :bin_min_frequency => 0.0 # Keep all bins for a single spectrum
-    )
-
-    # 3. Define the on_stage callback
-    collected_stage_data = []
-    stage_counter = Ref(0) # Initialize counter for sequential naming
-    normalized_spectrum_total = nothing # Variable to hold the normalized spectrum
-
-    function stage_callback_total(stage; idx, mz, intensity)
-        stage_counter[] += 1 # Increment counter
-        println("  -> Generating plot for stage: $stage (Total Spectrum)")
-        
-        local fig
-
-        if stage == :normalize
-            normalized_spectrum_total = (mz, intensity)
-            fig = plot_stage_spectrum(mz, intensity, title="Stage: $stage (Total Spectrum)")
-        elseif stage == :peaks && normalized_spectrum_total !== nothing
-            fig = Figure(size = (1400, 500))
-            ax = Axis(fig[1, 1], title="Stage: Peaks (Total Spectrum)", xlabel="m/z", ylabel="Intensity")
-            lines!(ax, normalized_spectrum_total[1], normalized_spectrum_total[2], color=:gray, label="Normalized Spectrum")
-            scatter!(ax, mz, intensity, color=:red, marker=:circle, markersize=8, label="Detected Peaks")
-            axislegend(ax)
-        else
-            fig = plot_stage_spectrum(mz, intensity, title="Stage: $stage (Total Spectrum)")
-        end
-
-        # Save separate plot with sequential name
-        stage_output_path = joinpath(output_dir, "$(file_type_prefix)_total_$(stage_counter[])_$(stage).png")
-        save(stage_output_path, fig)
-
-        # Collect data for overlay plot
-        push!(collected_stage_data, (stage, mz, intensity))
-    end
-
-    # 4. Run the pipeline on the single total spectrum
-    println("Running pipeline with steps: $pipeline_steps")
-    processed_result = run_preprocessing_pipeline(
-        [total_spectrum], # Pass the total spectrum as a vector of one spectrum
-        steps=pipeline_steps,
-        params=params,
-        on_stage=stage_callback_total
-    )
-
-    # 5. Generate and save the overlay plot
-    overlay_output_path = joinpath(output_dir, "$(file_type_prefix)_total_all_stages_overlay.png")
-    plot_overlay_stages(collected_stage_data, overlay_output_path, "Preprocessing Stages Overlay (Total Spectrum)")
-
-    # 6. Save feature matrix if generated
-    if processed_result isa FeatureMatrix
-        feature_matrix_output_path = joinpath(output_dir, "$(file_type_prefix)_total_feature_matrix.csv")
-        # Convert mz_bins to a more readable format for CSV
-        mz_labels = ["$(round(b[1], digits=4))_$(round(b[2], digits=4))" for b in processed_result.mz_bins]
-        df = DataFrame(processed_result.matrix, Symbol.(mz_labels))
-        CSV.write(feature_matrix_output_path, df)
-        println("SUCCESS: Feature matrix saved to $feature_matrix_output_path")
-    else
-        @warn "Pipeline did not return a FeatureMatrix for Total Spectrum."
-        processed_result
-    end
-
-    println("--- Pipeline test finished for TOTAL Spectrum (File Type: $file_type_prefix) ---")
-    println("Check the '$(output_dir)' directory for output plots and CSVs.")
-end
-
-
-# ===================================================================
-# TEST RUNNER
-# ===================================================================
-
-function run_preprocessing_tests()
-    println("="^80)
-    println("STARTING PREPROCESSING TEST SUITE")
-    println("="^80)
-
-    # --- Test Case 1: Run full pipeline on a single mzML spectrum ---
-    println("\n" * "="^20 * " Test Case 1: Full Pipeline on .mzML Spectrum " * "="^20)
-    println("FILE: ", TEST_MZML_FILE)
-    if isfile(TEST_MZML_FILE)
-        try
-            msi_data_mzml = OpenMSIData(TEST_MZML_FILE)
-            
-            # Dynamically determine tolerance
-            println("\n--- Calculating optimal tolerance for .mzML data ---")
-            report_mzml = analyze_mass_accuracy(msi_data_mzml, get_common_calibration_standards(:maldi_pos))
-            mz_tolerance_mzml = 0.002 # Default
-            if haskey(report_mzml, :optimal_ppm) && !isnan(report_mzml.optimal_ppm) && !isempty(report_mzml.matched_peaks)
-                avg_mz = mean([p[1] for p in report_mzml.matched_peaks])
-                mz_tolerance_mzml = avg_mz * report_mzml.optimal_ppm / 1e6
-                println("Optimal PPM: $(round(report_mzml.optimal_ppm, digits=2)), Average m/z: $(round(avg_mz, digits=2))")
-                println("Calculated m/z tolerance: $(round(mz_tolerance_mzml, digits=5))")
-            else
-                println("Could not determine optimal tolerance, using default: $mz_tolerance_mzml")
-            end
-
-            # Create a dedicated subdirectory for the output plots
-            mzml_output_dir = joinpath(RESULTS_DIR, "mzml_pipeline_stages")
-            mkpath(mzml_output_dir)
-            
-            test_full_pipeline(msi_data_mzml, MZML_SPECTRUM_ID, output_dir=mzml_output_dir, file_type_prefix="mzml", mz_tolerance=mz_tolerance_mzml)
-            test_full_pipeline_on_total_spectrum(msi_data_mzml, output_dir=mzml_output_dir, file_type_prefix="mzml", mz_tolerance=mz_tolerance_mzml)
-        catch e
-            println("ERROR in .mzML pipeline test: $e")
-            showerror(stdout, e, catch_backtrace())
-        end
-    else
-        println("SKIPPED: File not found: $TEST_MZML_FILE")
-    end
-
-    # --- Test Case 2: Run full pipeline on a single imzML spectrum ---
-    println("\n" * "="^20 * " Test Case 2: Full Pipeline on .imzML Spectrum " * "="^20)
-    println("FILE: ", TEST_IMZML_FILE)
-    if isfile(TEST_IMZML_FILE)
-        try
-            msi_data_imzml = OpenMSIData(TEST_IMZML_FILE)
-
-            # Dynamically determine tolerance
-            println("\n--- Calculating optimal tolerance for .imzML data ---")
-            report_imzml = analyze_mass_accuracy(msi_data_imzml, get_common_calibration_standards(:maldi_pos))
-            mz_tolerance_imzml = 0.002 # Default
-            if haskey(report_imzml, :optimal_ppm) && !isnan(report_imzml.optimal_ppm) && !isempty(report_imzml.matched_peaks)
-                avg_mz = mean([p[1] for p in report_imzml.matched_peaks])
-                mz_tolerance_imzml = avg_mz * report_imzml.optimal_ppm / 1e6
-                println("Optimal PPM: $(round(report_imzml.optimal_ppm, digits=2)), Average m/z: $(round(avg_mz, digits=2))")
-                println("Calculated m/z tolerance: $(round(mz_tolerance_imzml, digits=5))")
-            else
-                println("Could not determine optimal tolerance, using default: $mz_tolerance_imzml")
-            end
-            
-            # Create a dedicated subdirectory for the output plots
-            imzml_output_dir = joinpath(RESULTS_DIR, "imzml_pipeline_stages")
-            mkpath(imzml_output_dir)
-            
-            test_full_pipeline(msi_data_imzml, IMZML_COORDS, output_dir=imzml_output_dir, file_type_prefix="imzml", mz_tolerance=mz_tolerance_imzml)
-            test_full_pipeline_on_total_spectrum(msi_data_imzml, output_dir=imzml_output_dir, file_type_prefix="imzml", mz_tolerance=mz_tolerance_imzml)
-            # generate_qc_report(msi_data_imzml, TEST_IMZML_FILE, output_dir=imzml_output_dir)
-            custom_reference_peaks = Dict(
-                31.974 => "Red Phosphorus",
-                432.6584 => "P13",
-                464.6059 => "P15",
-                526.5534 => "P17",
-                650.4485 => "P21",
-                774.3435 => "P25",
-                898.2385 => "P29",
-                950.1861 => "P31",
-                1022.1336 => "P33",
-                1146.0286 => "P37",
-                1593.8187 => "P45",
-                772.433 => "Unknown 1",
-                772.5253 => "Unknown 2"
-            )
-            n_samples = length(msi_data_imzml.spectra_metadata)
-            generate_qc_report(msi_data_imzml, TEST_IMZML_FILE, reference_peaks=custom_reference_peaks, output_dir=imzml_output_dir, sample_spectra=n_samples)
-        catch e
-            println("ERROR in .imzML pipeline test: $e")
-            showerror(stdout, e, catch_backtrace())
-        end
-    else
-        println("SKIPPED: File not found: $TEST_IMZML_FILE")
-    end
-
-    # --- Test Case 3: Generate QC Report for .imzML data ---
-    println("\n" * "="^20 * " Test Case 3: QC Report Generation for .imzML " * "="^20)
-    println("FILE: ", TEST_IMZML_FILE)
-    if isfile(TEST_IMZML_FILE)
-        try
-            msi_data_imzml = OpenMSIData(TEST_IMZML_FILE)
-            
-            # Create a dedicated subdirectory for the QC report
-            qc_output_dir = joinpath(RESULTS_DIR, "qc_report")
-            mkpath(qc_output_dir)
-
-            println("\n--- Generating comprehensive QC report ---")
-            custom_reference_peaks = Dict(
-                31.974 => "Red Phosphorus",
-                432.6584 => "P13",
-                464.6059 => "P15",
-                526.5534 => "P17",
-                650.4485 => "P21",
-                774.3435 => "P25",
-                898.2385 => "P29",
-                950.1861 => "P31",
-                1022.1336 => "P33",
-                1146.0286 => "P37",
-                1593.8187 => "P45",
-                772.433 => "Unknown 1",
-                772.5253 => "Unknown 2"
-            )
-            # You can control the number of spectra sampled for the QC report.
-            # For the most accurate results, you can sample all spectra, but it will take longer.
-            # To sample all, use: n_samples = length(msi_data_imzml.spectra_metadata)
-            n_samples = length(msi_data_imzml.spectra_metadata)
-
-            #generate_qc_report(msi_data_imzml, TEST_IMZML_FILE, output_dir=qc_output_dir)
-            generate_qc_report(msi_data_imzml, TEST_IMZML_FILE, reference_peaks=custom_reference_peaks, output_dir=qc_output_dir, sample_spectra=n_samples)
-            
-            println("\n--- Analyzing specific reference peaks ---")
-            #=
-            reference_peaks = Dict(
-                104.10754 => "Imidazole",
-                175.11995 => "GPC fragment", 
-                226.15687 => "Phosphocholine"
-            )
-            =#
-            report = analyze_mass_accuracy(msi_data_imzml, custom_reference_peaks)
-            if haskey(report, :optimal_ppm)
-                println("Optimal PPM tolerance with specific peaks: $(round(report.optimal_ppm, digits=2)) ppm")
-            else
-                println("Could not determine optimal PPM with specific peaks.")
-            end
-
-        catch e
-            println("ERROR in QC report generation test: $e")
-            showerror(stdout, e, catch_backtrace())
-        end
-    else
-        println("SKIPPED: File not found: $TEST_IMZML_FILE")
-    end
-
-    println("\nPreprocessing tests finished.")
-end
 
 # --- Execute ---
-# Ensure the results directory exists
-mkpath(RESULTS_DIR)
-@time run_preprocessing_tests()
\ No newline at end of file
+run_preprocessing_suite()