diff --git a/Manifest.toml b/Manifest.toml
index 2547f94..8a1f9f5 100644
--- a/Manifest.toml
+++ b/Manifest.toml
@@ -1,14 +1,19 @@
 # This file is machine-generated - editing it directly is not advised
 
-julia_version = "1.7.3"
+julia_version = "1.9.3"
 manifest_format = "2.0"
-project_hash = "42cc740256e4fa7433ca5a6ab444ec66c69d7696"
+project_hash = "15f58ae716c5e6099f5c37ec8e385756d282b693"
 
 [[deps.AbstractFFTs]]
-deps = ["ChainRulesCore", "LinearAlgebra", "Test"]
+deps = ["LinearAlgebra"]
 git-tree-sha1 = "d92ad398961a3ed262d8bf04a1a2b8340f915fef"
 uuid = "621f4979-c628-5d54-868e-fcf4e3e8185c"
 version = "1.5.0"
+weakdeps = ["ChainRulesCore", "Test"]
+
+    [deps.AbstractFFTs.extensions]
+    AbstractFFTsChainRulesCoreExt = "ChainRulesCore"
+    AbstractFFTsTestExt = "Test"
 
 [[deps.AbstractLattices]]
 git-tree-sha1 = "f35684b7349da49fcc8a9e520e30e45dbb077166"
@@ -26,6 +31,18 @@ git-tree-sha1 = "954634616d5846d8e216df1298be2298d55280b2"
 uuid = "7d9f7c33-5ae7-4f3b-8dc6-eff91059b697"
 version = "0.1.32"
 
+    [deps.Accessors.extensions]
+    AccessorsAxisKeysExt = "AxisKeys"
+    AccessorsIntervalSetsExt = "IntervalSets"
+    AccessorsStaticArraysExt = "StaticArrays"
+    AccessorsStructArraysExt = "StructArrays"
+
+    [deps.Accessors.weakdeps]
+    AxisKeys = "94b1ba4f-4ee9-5380-92f1-94cde586c3c5"
+    IntervalSets = "8197267c-284f-5f27-9208-e0e47529a953"
+    StaticArrays = "90137ffa-7385-5640-81b9-e52037218182"
+    StructArrays = "09ab397b-f2b6-538f-b94a-2f83cf4a842a"
+
 [[deps.AccurateArithmetic]]
 deps = ["LinearAlgebra", "Random", "VectorizationBase"]
 git-tree-sha1 = "07af26e8d08c211ef85918f3e25d4c0990d20d70"
@@ -37,6 +54,10 @@ deps = ["LinearAlgebra", "Requires"]
 git-tree-sha1 = "76289dc51920fdc6e0013c872ba9551d54961c24"
 uuid = "79e6a3ab-5dfb-504d-930d-738a2a938a0e"
 version = "3.6.2"
+weakdeps = ["StaticArrays"]
+
+    [deps.Adapt.extensions]
+    AdaptStaticArraysExt = "StaticArrays"
 
 [[deps.Animations]]
 deps = ["Colors"]
@@ -51,6 +72,7 @@ version = "2.3.0"
 
 [[deps.ArgTools]]
 uuid = "0dad84c5-d112-42e6-8d28-ef12dabb789f"
+version = "1.1.1"
 
 [[deps.ArnoldiMethod]]
 deps = ["LinearAlgebra", "Random", "StaticArrays"]
@@ -76,6 +98,22 @@ git-tree-sha1 = "f83ec24f76d4c8f525099b2ac475fc098138ec31"
 uuid = "4fba245c-0d91-5ea0-9b3e-6abc04ee57a9"
 version = "7.4.11"
 
+    [deps.ArrayInterface.extensions]
+    ArrayInterfaceBandedMatricesExt = "BandedMatrices"
+    ArrayInterfaceBlockBandedMatricesExt = "BlockBandedMatrices"
+    ArrayInterfaceCUDAExt = "CUDA"
+    ArrayInterfaceGPUArraysCoreExt = "GPUArraysCore"
+    ArrayInterfaceStaticArraysCoreExt = "StaticArraysCore"
+    ArrayInterfaceTrackerExt = "Tracker"
+
+    [deps.ArrayInterface.weakdeps]
+    BandedMatrices = "aae01518-5342-5314-be14-df237901396f"
+    BlockBandedMatrices = "ffab5731-97b5-5995-9138-79e8c1846df0"
+    CUDA = "052768ef-5323-5732-b1bb-66c8b64840ba"
+    GPUArraysCore = "46192b85-c4d5-4398-a991-12ede77f4527"
+    StaticArraysCore = "1e83bf80-4336-4d27-bf5d-d5a4f845583c"
+    Tracker = "9f7883ad-71c0-57eb-9f7f-b5c9e6d3789c"
+
 [[deps.ArrayInterfaceCore]]
 deps = ["LinearAlgebra", "SnoopPrecompile", "SparseArrays", "SuiteSparse"]
 git-tree-sha1 = "e5f08b5689b1aad068e01751889f2f615c7db36d"
@@ -144,6 +182,20 @@ git-tree-sha1 = "e28912ce94077686443433c2800104b061a827ed"
 uuid = "198e06fe-97b7-11e9-32a5-e1d131e6ad66"
 version = "0.3.39"
 
+    [deps.BangBang.extensions]
+    BangBangChainRulesCoreExt = "ChainRulesCore"
+    BangBangDataFramesExt = "DataFrames"
+    BangBangStaticArraysExt = "StaticArrays"
+    BangBangStructArraysExt = "StructArrays"
+    BangBangTypedTablesExt = "TypedTables"
+
+    [deps.BangBang.weakdeps]
+    ChainRulesCore = "d360d2e6-b24c-11e9-a2a3-2a2ae2dbcce4"
+    DataFrames = "a93c6f00-e57d-5684-b7b6-d8193f3e46c0"
+    StaticArrays = "90137ffa-7385-5640-81b9-e52037218182"
+    StructArrays = "09ab397b-f2b6-538f-b94a-2f83cf4a842a"
+    TypedTables = "9d95f2ec-7b3d-5a63-8d20-e2491e220bb9"
+
 [[deps.Base64]]
 uuid = "2a0f44e3-6c83-55bd-87e4-b1978d98bd5f"
 
@@ -267,6 +319,18 @@ git-tree-sha1 = "1568b28f91293458345dabba6a5ea3f183250a61"
 uuid = "324d7699-5711-5eae-9e2f-1d82baa6b597"
 version = "0.10.8"
 
+    [deps.CategoricalArrays.extensions]
+    CategoricalArraysJSONExt = "JSON"
+    CategoricalArraysRecipesBaseExt = "RecipesBase"
+    CategoricalArraysSentinelArraysExt = "SentinelArrays"
+    CategoricalArraysStructTypesExt = "StructTypes"
+
+    [deps.CategoricalArrays.weakdeps]
+    JSON = "682c06a0-de6a-54ab-a142-c8b1cf79cde6"
+    RecipesBase = "3cdcf5f2-1ef4-517c-9805-6587b60abb01"
+    SentinelArrays = "91c51154-3ec4-41a3-a24f-3f23e20d615c"
+    StructTypes = "856f2bd8-1eba-4b0a-8007-ebc267875bd4"
+
 [[deps.Chain]]
 git-tree-sha1 = "8c4920235f6c561e401dfe569beb8b924adad003"
 uuid = "8be319e6-bccf-4806-a6f7-6fae938471bc"
@@ -279,10 +343,14 @@ uuid = "d360d2e6-b24c-11e9-a2a3-2a2ae2dbcce4"
 version = "1.16.0"
 
 [[deps.ChangesOfVariables]]
-deps = ["InverseFunctions", "LinearAlgebra", "Test"]
+deps = ["LinearAlgebra", "Test"]
 git-tree-sha1 = "2fba81a302a7be671aefe194f0525ef231104e7f"
 uuid = "9e997f8a-9a97-42d5-a9f1-ce6bfc15e2c0"
 version = "0.1.8"
+weakdeps = ["InverseFunctions"]
+
+    [deps.ChangesOfVariables.extensions]
+    ChangesOfVariablesInverseFunctionsExt = "InverseFunctions"
 
 [[deps.CloseOpenIntervals]]
 deps = ["Static", "StaticArrayInterface"]
@@ -373,19 +441,28 @@ uuid = "bbf7d656-a473-5ed7-a52c-81e309532950"
 version = "0.3.0"
 
 [[deps.Compat]]
-deps = ["Dates", "LinearAlgebra", "UUIDs"]
+deps = ["UUIDs"]
 git-tree-sha1 = "e460f044ca8b99be31d35fe54fc33a5c33dd8ed7"
 uuid = "34da2185-b29b-5c13-b0c7-acf172513d20"
 version = "4.9.0"
+weakdeps = ["Dates", "LinearAlgebra"]
+
+    [deps.Compat.extensions]
+    CompatLinearAlgebraExt = "LinearAlgebra"
 
 [[deps.CompilerSupportLibraries_jll]]
 deps = ["Artifacts", "Libdl"]
 uuid = "e66e0078-7015-5450-92f7-15fbd957f2ae"
+version = "1.0.5+0"
 
 [[deps.CompositionsBase]]
 git-tree-sha1 = "802bb88cd69dfd1509f6670416bd4434015693ad"
 uuid = "a33af91c-f02d-484b-be07-31d278c5ca2b"
 version = "0.1.2"
+weakdeps = ["InverseFunctions"]
+
+    [deps.CompositionsBase.extensions]
+    CompositionsBaseInverseFunctionsExt = "InverseFunctions"
 
 [[deps.ComputationalResources]]
 git-tree-sha1 = "52cb3ec90e8a8bea0e62e275ba577ad0f74821f7"
@@ -409,6 +486,11 @@ deps = ["LinearAlgebra"]
 git-tree-sha1 = "fe2838a593b5f776e1597e086dcd47560d94e816"
 uuid = "187b0558-2788-49d3-abe0-74a17ed4e7c9"
 version = "1.5.3"
+weakdeps = ["IntervalSets", "StaticArrays"]
+
+    [deps.ConstructionBase.extensions]
+    ConstructionBaseIntervalSetsExt = "IntervalSets"
+    ConstructionBaseStaticArraysExt = "StaticArrays"
 
 [[deps.Contour]]
 git-tree-sha1 = "d05d9e7b7aedff4e5b51a029dced05cfb6125781"
@@ -481,14 +563,20 @@ uuid = "244e2a9f-e319-4986-a169-4d1fe445cd52"
 version = "0.1.2"
 
 [[deps.DelaunayTriangulation]]
-deps = ["DataStructures", "EnumX", "ExactPredicates", "MakieCore", "Random", "SimpleGraphs"]
+deps = ["DataStructures", "EnumX", "ExactPredicates", "Random", "SimpleGraphs"]
 git-tree-sha1 = "a1d8532de83f8ce964235eff1edeff9581144d02"
 uuid = "927a84f5-c5f4-47a5-9785-b46e178433df"
 version = "0.7.2"
+weakdeps = ["MakieCore"]
+
+    [deps.DelaunayTriangulation.extensions]
+    DelaunayTriangulationMakieCoreExt = "MakieCore"
 
 [[deps.DelimitedFiles]]
 deps = ["Mmap"]
+git-tree-sha1 = "9e2f36d3c96a820c678f2f1f1782582fcf685bae"
 uuid = "8bb1440f-4735-579b-a4ab-409b98df4dab"
+version = "1.9.1"
 
 [[deps.DensityInterface]]
 deps = ["InverseFunctions", "Test"]
@@ -509,20 +597,29 @@ uuid = "b552c78f-8df3-52c6-915a-8e097449b14b"
 version = "1.15.1"
 
 [[deps.Distances]]
-deps = ["LinearAlgebra", "SparseArrays", "Statistics", "StatsAPI"]
+deps = ["LinearAlgebra", "Statistics", "StatsAPI"]
 git-tree-sha1 = "b6def76ffad15143924a2199f72a5cd883a2e8a9"
 uuid = "b4f34e82-e78d-54a5-968a-f98e89d6e8f7"
 version = "0.10.9"
+weakdeps = ["SparseArrays"]
+
+    [deps.Distances.extensions]
+    DistancesSparseArraysExt = "SparseArrays"
 
 [[deps.Distributed]]
 deps = ["Random", "Serialization", "Sockets"]
 uuid = "8ba89e20-285c-5b6f-9357-94700520ee1b"
 
 [[deps.Distributions]]
-deps = ["ChainRulesCore", "DensityInterface", "FillArrays", "LinearAlgebra", "PDMats", "Printf", "QuadGK", "Random", "SparseArrays", "SpecialFunctions", "Statistics", "StatsAPI", "StatsBase", "StatsFuns", "Test"]
+deps = ["FillArrays", "LinearAlgebra", "PDMats", "Printf", "QuadGK", "Random", "SparseArrays", "SpecialFunctions", "Statistics", "StatsAPI", "StatsBase", "StatsFuns", "Test"]
 git-tree-sha1 = "e76a3281de2719d7c81ed62c6ea7057380c87b1d"
 uuid = "31c24e10-a181-5473-b8eb-7969acd0382f"
 version = "0.25.98"
+weakdeps = ["ChainRulesCore", "DensityInterface"]
+
+    [deps.Distributions.extensions]
+    DistributionsChainRulesCoreExt = "ChainRulesCore"
+    DistributionsDensityInterfaceExt = "DensityInterface"
 
 [[deps.DocStringExtensions]]
 deps = ["LibGit2"]
@@ -539,6 +636,7 @@ version = "1.2.4"
 [[deps.Downloads]]
 deps = ["ArgTools", "FileWatching", "LibCURL", "NetworkOptions"]
 uuid = "f43a241f-c20a-4ad4-852c-f6b1247861c6"
+version = "1.6.0"
 
 [[deps.DualNumbers]]
 deps = ["Calculus", "NaNMath", "SpecialFunctions"]
@@ -557,6 +655,12 @@ git-tree-sha1 = "bdb1942cd4c45e3c678fd11569d5cccd80976237"
 uuid = "4e289a0a-7415-4d19-859d-a7e5c4648b56"
 version = "1.0.4"
 
+[[deps.EpollShim_jll]]
+deps = ["Artifacts", "JLLWrappers", "Libdl"]
+git-tree-sha1 = "8e9441ee83492030ace98f9789a654a6d0b1f643"
+uuid = "2702e6a9-849d-5ed8-8c21-79e8b8f9ee43"
+version = "0.0.20230411+0"
+
 [[deps.ErrorfreeArithmetic]]
 git-tree-sha1 = "d6863c556f1142a061532e79f611aa46be201686"
 uuid = "90fa49ef-747e-5e6f-a989-263ba693cf1a"
@@ -653,11 +757,21 @@ uuid = "1a297f60-69ca-5386-bcde-b61e274b549b"
 version = "1.5.0"
 
 [[deps.FiniteDiff]]
-deps = ["ArrayInterface", "LinearAlgebra", "Requires", "Setfield", "SparseArrays", "StaticArrays"]
+deps = ["ArrayInterface", "LinearAlgebra", "Requires", "Setfield", "SparseArrays"]
 git-tree-sha1 = "c6e4a1fbe73b31a3dea94b1da449503b8830c306"
 uuid = "6a86dc24-6348-571c-b903-95158fe2bd41"
 version = "2.21.1"
 
+    [deps.FiniteDiff.extensions]
+    FiniteDiffBandedMatricesExt = "BandedMatrices"
+    FiniteDiffBlockBandedMatricesExt = "BlockBandedMatrices"
+    FiniteDiffStaticArraysExt = "StaticArrays"
+
+    [deps.FiniteDiff.weakdeps]
+    BandedMatrices = "aae01518-5342-5314-be14-df237901396f"
+    BlockBandedMatrices = "ffab5731-97b5-5995-9138-79e8c1846df0"
+    StaticArrays = "90137ffa-7385-5640-81b9-e52037218182"
+
 [[deps.FixedPointNumbers]]
 deps = ["Statistics"]
 git-tree-sha1 = "335bfdceacc84c5cdf16aadc768aa5ddfc5383cc"
@@ -683,10 +797,14 @@ uuid = "59287772-0a20-5a39-b81b-1366585eb4c0"
 version = "0.4.2"
 
 [[deps.ForwardDiff]]
-deps = ["CommonSubexpressions", "DiffResults", "DiffRules", "LinearAlgebra", "LogExpFunctions", "NaNMath", "Preferences", "Printf", "Random", "SpecialFunctions", "StaticArrays"]
+deps = ["CommonSubexpressions", "DiffResults", "DiffRules", "LinearAlgebra", "LogExpFunctions", "NaNMath", "Preferences", "Printf", "Random", "SpecialFunctions"]
 git-tree-sha1 = "00e252f4d706b3d55a8863432e742bf5717b498d"
 uuid = "f6369f11-7733-5829-9624-2563aa707210"
 version = "0.10.35"
+weakdeps = ["StaticArrays"]
+
+    [deps.ForwardDiff.extensions]
+    ForwardDiffStaticArraysExt = "StaticArrays"
 
 [[deps.FreeType]]
 deps = ["CEnum", "FreeType2_jll"]
@@ -807,9 +925,9 @@ version = "1.3.14+0"
 
 [[deps.Graphs]]
 deps = ["ArnoldiMethod", "Compat", "DataStructures", "Distributed", "Inflate", "LinearAlgebra", "Random", "SharedArrays", "SimpleTraits", "SparseArrays", "Statistics"]
-git-tree-sha1 = "1cf1d7dcb4bc32d7b4a5add4232db3750c27ecb4"
+git-tree-sha1 = "899050ace26649433ef1af25bc17a815b3db52b7"
 uuid = "86223c79-3864-5bf0-83f7-82e725a168b6"
-version = "1.8.0"
+version = "1.9.0"
 
 [[deps.Grep]]
 deps = ["Test"]
@@ -1031,10 +1149,14 @@ uuid = "d1acc4aa-44c8-5952-acd4-ba5d80a2a253"
 version = "0.20.9"
 
 [[deps.IntervalSets]]
-deps = ["Dates", "Random", "Statistics"]
+deps = ["Dates", "Random"]
 git-tree-sha1 = "8e59ea773deee525c99a8018409f64f19fb719e6"
 uuid = "8197267c-284f-5f27-9208-e0e47529a953"
 version = "0.7.7"
+weakdeps = ["Statistics"]
+
+    [deps.IntervalSets.extensions]
+    IntervalSetsStatisticsExt = "Statistics"
 
 [[deps.InverseFunctions]]
 deps = ["Test"]
@@ -1168,6 +1290,14 @@ git-tree-sha1 = "f428ae552340899a935973270b8d98e5a31c49fe"
 uuid = "23fbe1c1-3f47-55db-b15f-69d7ec21a316"
 version = "0.16.1"
 
+    [deps.Latexify.extensions]
+    DataFramesExt = "DataFrames"
+    SymEngineExt = "SymEngine"
+
+    [deps.Latexify.weakdeps]
+    DataFrames = "a93c6f00-e57d-5684-b7b6-d8193f3e46c0"
+    SymEngine = "123dc426-2d89-5057-bbad-38513e3affd8"
+
 [[deps.LayoutPointers]]
 deps = ["ArrayInterface", "LinearAlgebra", "ManualMemory", "SIMDTypes", "Static", "StaticArrayInterface"]
 git-tree-sha1 = "88b8f66b604da079a627b6fb2860d3704a6729a1"
@@ -1200,10 +1330,12 @@ version = "0.1.0"
 [[deps.LibCURL]]
 deps = ["LibCURL_jll", "MozillaCACerts_jll"]
 uuid = "b27032c2-a3e7-50c8-80cd-2d36dbcbfd21"
+version = "0.6.3"
 
 [[deps.LibCURL_jll]]
 deps = ["Artifacts", "LibSSH2_jll", "Libdl", "MbedTLS_jll", "Zlib_jll", "nghttp2_jll"]
 uuid = "deac9b47-8bc7-5906-a0fe-35ac56dc84c0"
+version = "7.84.0+0"
 
 [[deps.LibGit2]]
 deps = ["Base64", "NetworkOptions", "Printf", "SHA"]
@@ -1212,6 +1344,7 @@ uuid = "76f85450-5226-5b5a-8eaa-529ad045b433"
 [[deps.LibSSH2_jll]]
 deps = ["Artifacts", "Libdl", "MbedTLS_jll"]
 uuid = "29816b5a-b9ab-546f-933c-edad1886dfa8"
+version = "1.10.2+0"
 
 [[deps.Libdl]]
 uuid = "8f399da3-3557-5675-b5ff-fb832c97cbdb"
@@ -1277,7 +1410,7 @@ uuid = "d3d80556-e9d4-5f37-9878-2ab0fcc64255"
 version = "7.2.0"
 
 [[deps.LinearAlgebra]]
-deps = ["Libdl", "libblastrampoline_jll"]
+deps = ["Libdl", "OpenBLAS_jll", "libblastrampoline_jll"]
 uuid = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 
 [[deps.LinearAlgebraX]]
@@ -1293,10 +1426,16 @@ uuid = "4345ca2d-374a-55d4-8d30-97f9976e7612"
 version = "0.6.1"
 
 [[deps.LogExpFunctions]]
-deps = ["ChainRulesCore", "ChangesOfVariables", "DocStringExtensions", "InverseFunctions", "IrrationalConstants", "LinearAlgebra"]
+deps = ["DocStringExtensions", "IrrationalConstants", "LinearAlgebra"]
 git-tree-sha1 = "c3ce8e7420b3a6e071e0fe4745f5d4300e37b13f"
 uuid = "2ab3a3ac-af41-5b50-aa03-7779005ae688"
 version = "0.3.24"
+weakdeps = ["ChainRulesCore", "ChangesOfVariables", "InverseFunctions"]
+
+    [deps.LogExpFunctions.extensions]
+    LogExpFunctionsChainRulesCoreExt = "ChainRulesCore"
+    LogExpFunctionsChangesOfVariablesExt = "ChangesOfVariables"
+    LogExpFunctionsInverseFunctionsExt = "InverseFunctions"
 
 [[deps.Logging]]
 uuid = "56ddb016-857b-54e1-b83d-db4d58db5568"
@@ -1308,10 +1447,15 @@ uuid = "e6f89c97-d47a-5376-807f-9c37f3926c36"
 version = "1.0.0"
 
 [[deps.LoopVectorization]]
-deps = ["ArrayInterface", "ArrayInterfaceCore", "CPUSummary", "ChainRulesCore", "CloseOpenIntervals", "DocStringExtensions", "ForwardDiff", "HostCPUFeatures", "IfElse", "LayoutPointers", "LinearAlgebra", "OffsetArrays", "PolyesterWeave", "PrecompileTools", "SIMDTypes", "SLEEFPirates", "SpecialFunctions", "Static", "StaticArrayInterface", "ThreadingUtilities", "UnPack", "VectorizationBase"]
+deps = ["ArrayInterface", "ArrayInterfaceCore", "CPUSummary", "CloseOpenIntervals", "DocStringExtensions", "HostCPUFeatures", "IfElse", "LayoutPointers", "LinearAlgebra", "OffsetArrays", "PolyesterWeave", "PrecompileTools", "SIMDTypes", "SLEEFPirates", "Static", "StaticArrayInterface", "ThreadingUtilities", "UnPack", "VectorizationBase"]
 git-tree-sha1 = "c88a4afe1703d731b1c4fdf4e3c7e77e3b176ea2"
 uuid = "bdcacae8-1622-11e9-2a5c-532679323890"
 version = "0.12.165"
+weakdeps = ["ChainRulesCore", "ForwardDiff", "SpecialFunctions"]
+
+    [deps.LoopVectorization.extensions]
+    ForwardDiffExt = ["ChainRulesCore", "ForwardDiff"]
+    SpecialFunctionsExt = "SpecialFunctions"
 
 [[deps.LsqFit]]
 deps = ["Distributions", "ForwardDiff", "LinearAlgebra", "NLSolversBase", "OptimBase", "Random", "StatsBase"]
@@ -1401,6 +1545,7 @@ version = "1.1.7"
 [[deps.MbedTLS_jll]]
 deps = ["Artifacts", "Libdl"]
 uuid = "c8ffd9c3-330d-5841-b78e-0817d7145fa1"
+version = "2.28.2+0"
 
 [[deps.Measures]]
 git-tree-sha1 = "c13304c81eec1ed3af7fc20e75fb6b26092a1102"
@@ -1441,6 +1586,7 @@ version = "0.3.4"
 
 [[deps.MozillaCACerts_jll]]
 uuid = "14a3606d-f60d-562e-9121-12d972cd8159"
+version = "2022.10.11"
 
 [[deps.MultipleTesting]]
 deps = ["Distributions", "SpecialFunctions", "StatsBase"]
@@ -1509,6 +1655,7 @@ version = "1.0.1"
 
 [[deps.NetworkOptions]]
 uuid = "ca575930-c2e3-43a9-ace4-1e988b2c1908"
+version = "1.2.0"
 
 [[deps.NodeJS]]
 deps = ["Pkg"]
@@ -1541,10 +1688,12 @@ version = "1.3.5+1"
 [[deps.OpenBLAS_jll]]
 deps = ["Artifacts", "CompilerSupportLibraries_jll", "Libdl"]
 uuid = "4536629a-c528-5b80-bd46-f80d51c5b363"
+version = "0.3.21+4"
 
 [[deps.OpenLibm_jll]]
 deps = ["Artifacts", "Libdl"]
 uuid = "05823500-19ac-5b8b-9628-191a04bc5112"
+version = "0.8.1+0"
 
 [[deps.OpenSSL_jll]]
 deps = ["Artifacts", "JLLWrappers", "Libdl"]
@@ -1584,6 +1733,7 @@ version = "1.6.2"
 [[deps.PCRE2_jll]]
 deps = ["Artifacts", "Libdl"]
 uuid = "efcefdf7-47ab-520b-bdef-62a2eaa19f15"
+version = "10.42.0+0"
 
 [[deps.PDMats]]
 deps = ["LinearAlgebra", "SparseArrays", "SuiteSparse"]
@@ -1651,8 +1801,9 @@ uuid = "30392449-352a-5448-841d-b1acce4e97dc"
 version = "0.42.2+0"
 
 [[deps.Pkg]]
-deps = ["Artifacts", "Dates", "Downloads", "LibGit2", "Libdl", "Logging", "Markdown", "Printf", "REPL", "Random", "SHA", "Serialization", "TOML", "Tar", "UUIDs", "p7zip_jll"]
+deps = ["Artifacts", "Dates", "Downloads", "FileWatching", "LibGit2", "Libdl", "Logging", "Markdown", "Printf", "REPL", "Random", "SHA", "Serialization", "TOML", "Tar", "UUIDs", "p7zip_jll"]
 uuid = "44cfe95a-1eb2-52ea-b672-e2afdf69b78f"
+version = "1.9.2"
 
 [[deps.PlotThemes]]
 deps = ["PlotUtils", "Statistics"]
@@ -1684,6 +1835,20 @@ git-tree-sha1 = "9f8675a55b37a70aa23177ec110f6e3f4dd68466"
 uuid = "91a5bcdd-55d7-5caf-9e0b-520d859cae80"
 version = "1.38.17"
 
+    [deps.Plots.extensions]
+    FileIOExt = "FileIO"
+    GeometryBasicsExt = "GeometryBasics"
+    IJuliaExt = "IJulia"
+    ImageInTerminalExt = "ImageInTerminal"
+    UnitfulExt = "Unitful"
+
+    [deps.Plots.weakdeps]
+    FileIO = "5789e2e9-d7fb-5bc7-8068-2c6fae9b9549"
+    GeometryBasics = "5c1252a2-5f33-56bf-86c9-59e7332b4326"
+    IJulia = "7073ff75-c697-5162-941a-fcdaad2a7d2a"
+    ImageInTerminal = "d8c32880-2388-543b-8c61-d9f865259254"
+    Unitful = "1986cc42-f94f-5a68-af5c-568840ba703d"
+
 [[deps.PolyesterWeave]]
 deps = ["BitTwiddlingConvenienceFunctions", "CPUSummary", "IfElse", "Static", "ThreadingUtilities"]
 git-tree-sha1 = "240d7170f5ffdb285f9427b92333c3463bf65bf6"
@@ -1696,10 +1861,16 @@ uuid = "647866c9-e3ac-4575-94e7-e3d426903924"
 version = "0.1.2"
 
 [[deps.Polynomials]]
-deps = ["ChainRulesCore", "LinearAlgebra", "MakieCore", "MutableArithmetics", "RecipesBase"]
+deps = ["LinearAlgebra", "RecipesBase"]
 git-tree-sha1 = "3aa2bb4982e575acd7583f01531f241af077b163"
 uuid = "f27b6e38-b328-58d1-80ce-0feddd5e7a45"
 version = "3.2.13"
+weakdeps = ["ChainRulesCore", "MakieCore", "MutableArithmetics"]
+
+    [deps.Polynomials.extensions]
+    PolynomialsChainRulesCoreExt = "ChainRulesCore"
+    PolynomialsMakieCoreExt = "MakieCore"
+    PolynomialsMutableArithmeticsExt = "MutableArithmetics"
 
 [[deps.PooledArrays]]
 deps = ["DataAPI", "Future"]
@@ -1807,6 +1978,10 @@ deps = ["Requires"]
 git-tree-sha1 = "1342a47bf3260ee108163042310d26f2be5ec90b"
 uuid = "c84ed2f1-dad5-54f0-aa8e-dbefe2724439"
 version = "0.4.5"
+weakdeps = ["FixedPointNumbers"]
+
+    [deps.Ratios.extensions]
+    RatiosFixedPointNumbersExt = "FixedPointNumbers"
 
 [[deps.RealDot]]
 deps = ["LinearAlgebra"]
@@ -1879,6 +2054,16 @@ git-tree-sha1 = "de432823e8aab4dd1a985be4be768f95acf152d4"
 uuid = "f2b01f46-fcfa-551c-844a-d8ac1e96c665"
 version = "2.0.17"
 
+    [deps.Roots.extensions]
+    RootsForwardDiffExt = "ForwardDiff"
+    RootsIntervalRootFindingExt = "IntervalRootFinding"
+    RootsSymPyExt = "SymPy"
+
+    [deps.Roots.weakdeps]
+    ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210"
+    IntervalRootFinding = "d2bf35a9-74e0-55ec-b149-d360ff49b807"
+    SymPy = "24249f21-da20-56a4-8eb1-6a02cf4ae2e6"
+
 [[deps.Rotations]]
 deps = ["LinearAlgebra", "Quaternions", "Random", "StaticArrays"]
 git-tree-sha1 = "54ccb4dbab4b1f69beb255a2c0ca5f65a9c82f08"
@@ -1892,6 +2077,7 @@ version = "0.2.1"
 
 [[deps.SHA]]
 uuid = "ea8e919c-243c-51af-8825-aaa63cd721ce"
+version = "0.7.0"
 
 [[deps.SIMD]]
 deps = ["PrecompileTools"]
@@ -1999,6 +2185,12 @@ git-tree-sha1 = "5d7e3f4e11935503d3ecaf7186eac40602e7d231"
 uuid = "699a6c99-e7fa-54fc-8d76-47d257e15c1d"
 version = "0.9.4"
 
+[[deps.SimpleWeightedGraphs]]
+deps = ["Graphs", "LinearAlgebra", "Markdown", "SparseArrays"]
+git-tree-sha1 = "4b33e0e081a825dbfaf314decf58fa47e53d6acb"
+uuid = "47aef6b3-ad0c-573a-a1e2-d07658019622"
+version = "1.4.0"
+
 [[deps.SnoopPrecompile]]
 deps = ["Preferences"]
 git-tree-sha1 = "e760a70afdcd461cf01a575947738d359234665c"
@@ -2027,14 +2219,18 @@ uuid = "a2af1166-a08f-5f64-846c-94a0d3cef48c"
 version = "1.1.1"
 
 [[deps.SparseArrays]]
-deps = ["LinearAlgebra", "Random"]
+deps = ["Libdl", "LinearAlgebra", "Random", "Serialization", "SuiteSparse_jll"]
 uuid = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"
 
 [[deps.SpecialFunctions]]
-deps = ["ChainRulesCore", "IrrationalConstants", "LogExpFunctions", "OpenLibm_jll", "OpenSpecFun_jll"]
+deps = ["IrrationalConstants", "LogExpFunctions", "OpenLibm_jll", "OpenSpecFun_jll"]
 git-tree-sha1 = "7beb031cf8145577fbccacd94b8a8f4ce78428d3"
 uuid = "276daf66-3868-5448-9aa4-cd146d93841b"
 version = "2.3.0"
+weakdeps = ["ChainRulesCore"]
+
+    [deps.SpecialFunctions.extensions]
+    SpecialFunctionsChainRulesCoreExt = "ChainRulesCore"
 
 [[deps.SplittablesBase]]
 deps = ["Setfield", "Test"]
@@ -2065,12 +2261,21 @@ deps = ["ArrayInterface", "Compat", "IfElse", "LinearAlgebra", "Requires", "Snoo
 git-tree-sha1 = "33040351d2403b84afce74dae2e22d3f5b18edcb"
 uuid = "0d7ed370-da01-4f52-bd93-41d350b8b718"
 version = "1.4.0"
+weakdeps = ["OffsetArrays", "StaticArrays"]
+
+    [deps.StaticArrayInterface.extensions]
+    StaticArrayInterfaceOffsetArraysExt = "OffsetArrays"
+    StaticArrayInterfaceStaticArraysExt = "StaticArrays"
 
 [[deps.StaticArrays]]
-deps = ["LinearAlgebra", "Random", "StaticArraysCore", "Statistics"]
+deps = ["LinearAlgebra", "Random", "StaticArraysCore"]
 git-tree-sha1 = "9cabadf6e7cd2349b6cf49f1915ad2028d65e881"
 uuid = "90137ffa-7385-5640-81b9-e52037218182"
 version = "1.6.2"
+weakdeps = ["Statistics"]
+
+    [deps.StaticArrays.extensions]
+    StaticArraysStatisticsExt = "Statistics"
 
 [[deps.StaticArraysCore]]
 git-tree-sha1 = "36b3d696ce6366023a0ea192b4cd442268995a0d"
@@ -2080,6 +2285,7 @@ version = "1.4.2"
 [[deps.Statistics]]
 deps = ["LinearAlgebra", "SparseArrays"]
 uuid = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"
+version = "1.9.0"
 
 [[deps.StatsAPI]]
 deps = ["LinearAlgebra"]
@@ -2094,10 +2300,15 @@ uuid = "2913bbd2-ae8a-5f71-8c99-4fb6c76f3a91"
 version = "0.33.21"
 
 [[deps.StatsFuns]]
-deps = ["ChainRulesCore", "HypergeometricFunctions", "InverseFunctions", "IrrationalConstants", "LogExpFunctions", "Reexport", "Rmath", "SpecialFunctions"]
+deps = ["HypergeometricFunctions", "IrrationalConstants", "LogExpFunctions", "Reexport", "Rmath", "SpecialFunctions"]
 git-tree-sha1 = "f625d686d5a88bcd2b15cd81f18f98186fdc0c9a"
 uuid = "4c63d2b9-4356-54db-8cca-17b64c39e42c"
 version = "1.3.0"
+weakdeps = ["ChainRulesCore", "InverseFunctions"]
+
+    [deps.StatsFuns.extensions]
+    StatsFunsChainRulesCoreExt = "ChainRulesCore"
+    StatsFunsInverseFunctionsExt = "InverseFunctions"
 
 [[deps.StatsModels]]
 deps = ["DataAPI", "DataStructures", "LinearAlgebra", "Printf", "REPL", "ShiftedArrays", "SparseArrays", "StatsBase", "StatsFuns", "Tables"]
@@ -2132,9 +2343,15 @@ version = "0.6.15"
 deps = ["Libdl", "LinearAlgebra", "Serialization", "SparseArrays"]
 uuid = "4607b0f0-06f3-5cda-b6b1-a6196a1729e9"
 
+[[deps.SuiteSparse_jll]]
+deps = ["Artifacts", "Libdl", "Pkg", "libblastrampoline_jll"]
+uuid = "bea87d4a-7f5b-5778-9afe-8cc45184846c"
+version = "5.10.1+6"
+
 [[deps.TOML]]
 deps = ["Dates"]
 uuid = "fa267f1f-6049-4f14-aa54-33bafae1ed76"
+version = "1.0.3"
 
 [[deps.TSne]]
 deps = ["Distances", "LinearAlgebra", "Printf", "ProgressMeter", "Statistics"]
@@ -2175,6 +2392,7 @@ version = "1.10.1"
 [[deps.Tar]]
 deps = ["ArgTools", "SHA"]
 uuid = "a4e569a6-e804-4fa4-b0f3-eef7a1d5b13e"
+version = "1.10.0"
 
 [[deps.TensorCore]]
 deps = ["LinearAlgebra"]
@@ -2228,6 +2446,20 @@ git-tree-sha1 = "53bd5978b182fa7c57577bdb452c35e5b4fb73a5"
 uuid = "28d57a85-8fef-5791-bfe6-a80928e7c999"
 version = "0.4.78"
 
+    [deps.Transducers.extensions]
+    TransducersBlockArraysExt = "BlockArrays"
+    TransducersDataFramesExt = "DataFrames"
+    TransducersLazyArraysExt = "LazyArrays"
+    TransducersOnlineStatsBaseExt = "OnlineStatsBase"
+    TransducersReferenceablesExt = "Referenceables"
+
+    [deps.Transducers.weakdeps]
+    BlockArrays = "8e7c35d0-a365-5155-bbbb-fb81a777f24e"
+    DataFrames = "a93c6f00-e57d-5684-b7b6-d8193f3e46c0"
+    LazyArrays = "5078a376-72f3-5289-bfd5-ec5146d43c02"
+    OnlineStatsBase = "925886fa-5bf2-5e8e-b522-a9147a512338"
+    Referenceables = "42d2dcc6-99eb-4e98-b66c-637b7d73030e"
+
 [[deps.TriplotBase]]
 git-tree-sha1 = "4d4ed7f294cda19382ff7de4c137d24d16adc89b"
 uuid = "981d1d27-644d-49a2-9326-4793e63143c3"
@@ -2274,10 +2506,15 @@ uuid = "1cfade01-22cf-5700-b092-accc4b62d6e1"
 version = "0.4.1"
 
 [[deps.Unitful]]
-deps = ["ConstructionBase", "Dates", "InverseFunctions", "LinearAlgebra", "Random"]
+deps = ["Dates", "LinearAlgebra", "Random"]
 git-tree-sha1 = "1cd9b6d3f637988ca788007b7466c132feebe263"
 uuid = "1986cc42-f94f-5a68-af5c-568840ba703d"
 version = "1.16.1"
+weakdeps = ["ConstructionBase", "InverseFunctions"]
+
+    [deps.Unitful.extensions]
+    ConstructionBaseUnitfulExt = "ConstructionBase"
+    InverseFunctionsUnitfulExt = "InverseFunctions"
 
 [[deps.UnitfulLatexify]]
 deps = ["LaTeXStrings", "Latexify", "Unitful"]
@@ -2320,10 +2557,14 @@ uuid = "81def892-9a0e-5fdd-b105-ffc91e053289"
 version = "1.3.0"
 
 [[deps.Wayland_jll]]
-deps = ["Artifacts", "Expat_jll", "JLLWrappers", "Libdl", "Libffi_jll", "Pkg", "XML2_jll"]
+deps = ["Artifacts", "EpollShim_jll", "Expat_jll", "JLLWrappers", "Libdl", "Libffi_jll", "Pkg", "XML2_jll"]
+<<<<<<< HEAD
 git-tree-sha1 = "ed8d92d9774b077c53e1da50fd81a36af3744c1c"
+=======
+git-tree-sha1 = "7558e29847e99bc3f04d6569e82d0f5c54460703"
+>>>>>>> 3d1a648 (updates)
 uuid = "a2964d1f-97da-50d4-b82a-358c7fce9d89"
-version = "1.21.0+0"
+version = "1.21.0+1"
 
 [[deps.Wayland_protocols_jll]]
 deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
@@ -2536,6 +2777,7 @@ version = "4.3.4+0"
 [[deps.Zlib_jll]]
 deps = ["Libdl"]
 uuid = "83775a58-1f1d-513f-b197-d71354ab007a"
+version = "1.2.13+0"
 
 [[deps.Zstd_jll]]
 deps = ["Artifacts", "JLLWrappers", "Libdl"]
@@ -2568,8 +2810,9 @@ uuid = "0ac62f75-1d6f-5e53-bd7c-93b484bb37c0"
 version = "0.15.1+0"
 
 [[deps.libblastrampoline_jll]]
-deps = ["Artifacts", "Libdl", "OpenBLAS_jll"]
+deps = ["Artifacts", "Libdl"]
 uuid = "8e850b90-86db-534c-a0d3-1478176c7d93"
+version = "5.8.0+0"
 
 [[deps.libfdk_aac_jll]]
 deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
@@ -2598,10 +2841,12 @@ version = "1.3.7+1"
 [[deps.nghttp2_jll]]
 deps = ["Artifacts", "Libdl"]
 uuid = "8e850ede-7688-5339-a07c-302acd2aaf8d"
+version = "1.48.0+0"
 
 [[deps.p7zip_jll]]
 deps = ["Artifacts", "Libdl"]
 uuid = "3f19e933-33d8-53b3-aaab-bd5110c3b7a0"
+version = "17.4.0+0"
 
 [[deps.x264_jll]]
 deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
@@ -2617,6 +2862,6 @@ version = "3.5.0+0"
 
 [[deps.xkbcommon_jll]]
 deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg", "Wayland_jll", "Wayland_protocols_jll", "Xorg_libxcb_jll", "Xorg_xkeyboard_config_jll"]
-git-tree-sha1 = "9ebfc140cc56e8c2156a15ceac2f0302e327ac0a"
+git-tree-sha1 = "9c304562909ab2bab0262639bd4f444d7bc2be37"
 uuid = "d8fb68d0-12a3-5cfd-a85a-d49703b185fd"
-version = "1.4.1+0"
+version = "1.4.1+1"
diff --git a/Project.toml b/Project.toml
index 70a7591..ed7bff5 100644
--- a/Project.toml
+++ b/Project.toml
@@ -16,20 +16,24 @@ DataFrames = "a93c6f00-e57d-5684-b7b6-d8193f3e46c0"
 DataFramesMeta = "1313f7d8-7da2-5740-9ea0-a2ca25f37964"
 Distances = "b4f34e82-e78d-54a5-968a-f98e89d6e8f7"
 Distributed = "8ba89e20-285c-5b6f-9357-94700520ee1b"
+FFTW = "7a1cc6ca-52ef-59f5-83cd-3a7055c09341"
 FileIO = "5789e2e9-d7fb-5bc7-8068-2c6fae9b9549"
 Folds = "41a02a25-b8f0-4f67-bc48-60067656b558"
 GZip = "92fee26a-97fe-5a0c-ad85-20a5f3185b63"
 GeneticsMakie = "8ca62643-82d8-47b5-a233-a06d1654fb35"
 Glob = "c27321d9-0574-5035-807b-f59d2c89b15c"
+Graphs = "86223c79-3864-5bf0-83f7-82e725a168b6"
 Grep = "775960c6-9b90-5df0-b405-1e337feb71e5"
 HDF5 = "f67ccb44-e63f-5c2f-98bd-6dc0ccc4ba2f"
 HTTP = "cd3eb016-35fb-5094-929b-558a96fad6f3"
 HypothesisTests = "09f84164-cd44-5f33-b23f-e6b0d136a0d5"
 IJulia = "7073ff75-c697-5162-941a-fcdaad2a7d2a"
 ImageCore = "a09fc81d-aa75-5fe9-8630-4744c3626534"
+ImageMorphology = "787d08f9-d448-5407-9aad-5290dd7ab264"
 Images = "916415d5-f1e6-5110-898d-aaa5f9f070e0"
 JLD2 = "033835bb-8acc-5ee8-8aae-3f567f8a3819"
 JSON = "682c06a0-de6a-54ab-a142-c8b1cf79cde6"
+KernelDensity = "5ab0869b-81aa-558d-bb23-cbf5423bbe9b"
 Leiden = "e5cf005d-d64d-44d4-821e-2384bf4ace83"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 Loess = "4345ca2d-374a-55d4-8d30-97f9976e7612"
@@ -38,6 +42,7 @@ MatrixMarket = "4d4711f2-db25-561a-b6b3-d35e7d4047d3"
 MultipleTesting = "f8716d33-7c4a-5097-896f-ce0ecbd3ef6b"
 MultivariateStats = "6f286f6a-111f-5878-ab1e-185364afe411"
 NearestNeighborDescent = "dd2c4c9e-a32f-5b2f-b342-08c2f244fce8"
+NearestNeighbors = "b8a86587-4115-5ab1-83bc-aa920d37bbce"
 Pkg = "44cfe95a-1eb2-52ea-b672-e2afdf69b78f"
 PlotlyJS = "f0f68f2c-4968-5e81-91da-67840de0976a"
 Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80"
@@ -48,6 +53,7 @@ PyCallUtils = "02651160-e6c3-11e9-0e2a-551f75d9519e"
 RCall = "6f49c342-dc21-5d91-9882-a32aef131414"
 Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
 RollingFunctions = "b0e4dd01-7b14-53d8-9b45-175a3e362653"
+SimpleWeightedGraphs = "47aef6b3-ad0c-573a-a1e2-d07658019622"
 SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"
 Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"
 StatsBase = "2913bbd2-ae8a-5f71-8c99-4fb6c76f3a91"
diff --git a/src/CellScopes.jl b/src/CellScopes.jl
index 858a816..b75122b 100644
--- a/src/CellScopes.jl
+++ b/src/CellScopes.jl
@@ -54,6 +54,14 @@ using GeneticsMakie
 using Arrow
 using HDF5
 using PyCallUtils
+import FFTW
+import ImageMorphology
+import Graphs
+using NearestNeighbors: KDTree, knn
+using SimpleWeightedGraphs
+using StatsBase: countmap
+using KernelDensity
+using Graphs: src, dst
 
 include("scrna/objects.jl")
 include("spatial/sp_objects.jl")
@@ -70,4 +78,6 @@ include("spatial/sp_processing.jl")
 include("scatac/atac_processing.jl")
 include("scatac/atac_plots.jl")
 include("scatac/atac_utils.jl")
+include("spatial/baysor_boundary.jl")
+
 end
diff --git a/src/properties.jl b/src/properties.jl
index dca5d2b..561637f 100644
--- a/src/properties.jl
+++ b/src/properties.jl
@@ -56,7 +56,7 @@ function Base.show(io::IO, count::AbstractCount)
     [println("- ", string(i)) for i in fieldnames(typeof(count))]
 end
 
-function Base.show(io::IO, obj_type::Union{AbstractDimReduction, ClusteringObject, VariableGeneObject, UndefinedObject, SpaImputeObj, VisiumImgObject, SpaCountObj,SpaMetaObj})
+function Base.show(io::IO, obj_type::Union{AbstractDimReduction, ClusteringObject, VariableGeneObject, UndefinedObject, SpaImputeObj, VisiumImgObject, SpaCountObj,SpaMetaObj,SpaImageObj})
     println(io, string(typeof(obj_type)))
     println("All fields:")
     [println("- ", string(i)) for i in fieldnames(typeof(obj_type))]
@@ -76,10 +76,11 @@ matchtype(field) = @match field begin
     spMeta::SpaMetaObj => println("- Spatial metadata")
     spCoord::SpaCoordObj => println("- Spatial coordiantes")
     fragment::FragmentObject => println("- Fragment data")
+    imgObject::SpaImageObj => println("- Image data")
     uns::UndefinedObject => println("- Undefined slot")
 end
 
-function Base.show(io::IO, sp_obj::Union{ImagingSpatialObject, CartanaObject, VisiumObject, XeniumObject, MerfishObject, SlideseqObject, STARmapObject, seqFishObject})
+function Base.show(io::IO, sp_obj::Union{ImagingSpatialObject, CartanaObject, VisiumObject, XeniumObject, MerfishObject, SlideseqObject, STARmapObject, seqFishObject, StereoSeqObject})
     if isa(sp_obj, ImagingSpatialObject)
         println(io, "ImagingSpatialObject in CellScopes.jl")
     elseif isa(sp_obj, CartanaObject)
@@ -94,6 +95,8 @@ function Base.show(io::IO, sp_obj::Union{ImagingSpatialObject, CartanaObject, Vi
         println(io, "seqFishObject in CellScopes.jl")
     elseif isa(sp_obj, SlideseqObject)
         println(io, "SlideseqObject in CellScopes.jl")
+    elseif isa(sp_obj, StereoSeqObject)
+        println(io, "StereoSeqObject in CellScopes.jl")
     else
         println(io, "VisiumObject in CellScopes.jl")
     end
diff --git a/src/scrna/processing.jl b/src/scrna/processing.jl
index 3dac359..54ddbb1 100644
--- a/src/scrna/processing.jl
+++ b/src/scrna/processing.jl
@@ -12,7 +12,7 @@ function normalize_object(ct_obj::RawCountObject; scale_factor = 10000, norm_met
     return norm_obj
 end
 
-function normalize_object(sc_obj::Union{scRNAObject, VisiumObject,ImagingSpatialObject, CartanaObject, XeniumObject, MerfishObject, SlideseqObject, STARmapObject, seqFishObject}; scale_factor = 10000, norm_method = "logarithm", pseudocount = 1)
+function normalize_object(sc_obj::Union{scRNAObject, VisiumObject,ImagingSpatialObject, CartanaObject, XeniumObject, MerfishObject, SlideseqObject, STARmapObject, seqFishObject, StereoSeqObject}; scale_factor = 10000, norm_method = "logarithm", pseudocount = 1)
     norm_obj = normalize_object(sc_obj.rawCount; scale_factor = scale_factor, norm_method = norm_method, pseudocount = pseudocount)
     sc_obj.normCount = norm_obj
     if isa(sc_obj, Union{MerfishObject, CartanaObject, XeniumObject})
@@ -42,7 +42,7 @@ function scale_object(ct_obj::NormCountObject; features::Union{Vector{String}, N
     return scale_obj
 end
 
-function scale_object(sc_obj::Union{scRNAObject, VisiumObject, ImagingSpatialObject, CartanaObject, XeniumObject,MerfishObject, SlideseqObject, STARmapObject, seqFishObject}; features::Union{Vector{String}, Nothing}=nothing, scale_max = 10.0, do_scale::Bool = true, do_center::Bool = true)
+function scale_object(sc_obj::Union{scRNAObject, VisiumObject, ImagingSpatialObject, CartanaObject, XeniumObject,MerfishObject, SlideseqObject, STARmapObject, seqFishObject, StereoSeqObject}; features::Union{Vector{String}, Nothing}=nothing, scale_max = 10.0, do_scale::Bool = true, do_center::Bool = true)
     scale_obj = scale_object(sc_obj.normCount; features = features, scale_max=scale_max, do_scale=do_scale, do_center=do_center)
     sc_obj.scaleCount = scale_obj
     return sc_obj
@@ -72,14 +72,14 @@ function find_variable_genes(ct_mtx::RawCountObject; nFeatures::Int64 = 2000, sp
     return vst_data, Features
 end
 
-function find_variable_genes(sc_obj::Union{scRNAObject, VisiumObject, ImagingSpatialObject, CartanaObject, XeniumObject,MerfishObject, SlideseqObject, STARmapObject, seqFishObject}; nFeatures::Int64 = 2000, span::Float64 = 0.3)
+function find_variable_genes(sc_obj::Union{scRNAObject, VisiumObject, ImagingSpatialObject, CartanaObject, XeniumObject,MerfishObject, SlideseqObject, STARmapObject, seqFishObject, StereoSeqObject}; nFeatures::Int64 = 2000, span::Float64 = 0.3)
     vst_data, Features = find_variable_genes(sc_obj.rawCount;  nFeatures = nFeatures, span = span)
     var_obj = VariableGeneObject(Features, vst_data)
     sc_obj.varGene = var_obj
     return sc_obj
 end
 
-function run_pca(sc_obj::Union{scRNAObject, VisiumObject, ImagingSpatialObject, CartanaObject, XeniumObject,MerfishObject, SlideseqObject, STARmapObject, seqFishObject}; method=:svd, pratio = 1, maxoutdim = 10)
+function run_pca(sc_obj::Union{scRNAObject, VisiumObject, ImagingSpatialObject, CartanaObject, XeniumObject,MerfishObject, SlideseqObject, STARmapObject, seqFishObject, StereoSeqObject}; method=:svd, pratio = 1, maxoutdim = 10)
     features = sc_obj.varGene.var_gene
     if length(sc_obj.scaleCount.gene_name) == length(sc_obj.rawCount.gene_name)
         new_count = subset_count(sc_obj.scaleCount; genes = features)
@@ -99,7 +99,7 @@ function run_pca(sc_obj::Union{scRNAObject, VisiumObject, ImagingSpatialObject,
     return sc_obj
 end
 
-function run_clustering_atlas(sc_obj::Union{scRNAObject, VisiumObject, ImagingSpatialObject, CartanaObject, XeniumObject, scATACObject,MerfishObject, SlideseqObject, STARmapObject, seqFishObject}; n_neighbors=30, metric=CosineDist(), res= 0.06, seed_use=1234)
+function run_clustering_atlas(sc_obj::Union{scRNAObject, VisiumObject, ImagingSpatialObject, CartanaObject, XeniumObject, scATACObject,MerfishObject, SlideseqObject, STARmapObject, seqFishObject, StereoSeqObject}; n_neighbors=30, metric=CosineDist(), res= 0.06, seed_use=1234)
     if isdefined(sc_obj.dimReduction, :umap)
         indices = sc_obj.dimReduction.umap.knn_data
     else
@@ -147,7 +147,7 @@ function run_clustering_atlas(sc_obj::Union{scRNAObject, VisiumObject, ImagingSp
     return sc_obj
 end
 
-function run_clustering_small(sc_obj::Union{scRNAObject, VisiumObject, ImagingSpatialObject, CartanaObject, XeniumObject, scATACObject, MerfishObject, SlideseqObject, STARmapObject, seqFishObject}; n_neighbors=30, metric=CosineDist(), res= 0.06, seed_use=1234)
+function run_clustering_small(sc_obj::Union{scRNAObject, VisiumObject, ImagingSpatialObject, CartanaObject, XeniumObject, scATACObject, MerfishObject, SlideseqObject, STARmapObject, seqFishObject, StereoSeqObject}; n_neighbors=30, metric=CosineDist(), res= 0.06, seed_use=1234)
     knn_data = [collect(i) for i in eachrow(sc_obj.dimReduction.pca.cell_embedding)]
     graph = nndescent(knn_data, n_neighbors, metric)
     indices, dist_mat = knn_matrices(graph);
@@ -192,7 +192,7 @@ function run_clustering_small(sc_obj::Union{scRNAObject, VisiumObject, ImagingSp
     return sc_obj
 end
 
-function run_clustering(sc_obj::Union{scRNAObject, VisiumObject, ImagingSpatialObject, CartanaObject, XeniumObject, scATACObject, MerfishObject, SlideseqObject, STARmapObject, seqFishObject}; n_neighbors=30, metric=CosineDist(), res= 0.06, seed_use=1234)
+function run_clustering(sc_obj::Union{scRNAObject, VisiumObject, ImagingSpatialObject, CartanaObject, XeniumObject, scATACObject, MerfishObject, SlideseqObject, STARmapObject, seqFishObject, StereoSeqObject}; n_neighbors=30, metric=CosineDist(), res= 0.06, seed_use=1234)
     n = size(sc_obj.rawCount.count_mtx, 2)
     if n < 10000
         obj = run_clustering_small(sc_obj; n_neighbors=n_neighbors, metric=metric, res= res, seed_use=seed_use)
@@ -203,7 +203,7 @@ function run_clustering(sc_obj::Union{scRNAObject, VisiumObject, ImagingSpatialO
     end
 end
 
-function run_tsne(sc_obj::Union{scRNAObject, VisiumObject, ImagingSpatialObject, CartanaObject, XeniumObject, scATACObject, MerfishObject, SlideseqObject, STARmapObject, seqFishObject}; ndim::Int64 = 2, dims_use = 1:10, max_iter::Int64 = 2000, perplexit::Real = 30.0, pca_init::Bool = true,  seed_use::Int64 = 1234)
+function run_tsne(sc_obj::Union{scRNAObject, VisiumObject, ImagingSpatialObject, CartanaObject, XeniumObject, scATACObject, MerfishObject, SlideseqObject, STARmapObject, seqFishObject, StereoSeqObject}; ndim::Int64 = 2, dims_use = 1:10, max_iter::Int64 = 2000, perplexit::Real = 30.0, pca_init::Bool = true,  seed_use::Int64 = 1234)
     Random.seed!(seed_use)
     pca_mat = sc_obj.dimReduction.pca.cell_embedding
     pca_mat = pca_mat[:, dims_use]
@@ -215,7 +215,7 @@ function run_tsne(sc_obj::Union{scRNAObject, VisiumObject, ImagingSpatialObject,
     return sc_obj
 end
 
-function run_umap(sc_obj::Union{scRNAObject, VisiumObject, ImagingSpatialObject, CartanaObject, XeniumObject, scATACObject,MerfishObject, SlideseqObject, STARmapObject, seqFishObject}; ndim::Int64 = 2, dims_use = 1:10, n_neighbors::Int64 = 30, n_epochs=300, init = :spectral, metric = CosineDist(), min_dist::Real = 0.4, seed_use::Int64 = 1234)
+function run_umap(sc_obj::Union{scRNAObject, VisiumObject, ImagingSpatialObject, CartanaObject, XeniumObject, scATACObject,MerfishObject, SlideseqObject, STARmapObject, seqFishObject, StereoSeqObject}; ndim::Int64 = 2, dims_use = 1:10, n_neighbors::Int64 = 30, n_epochs=300, init = :spectral, metric = CosineDist(), min_dist::Real = 0.4, seed_use::Int64 = 1234)
     Random.seed!(seed_use)
     pca_mat = sc_obj.dimReduction.pca.cell_embedding
     pca_mat = pca_mat[:, dims_use]
@@ -233,7 +233,7 @@ function run_umap(sc_obj::Union{scRNAObject, VisiumObject, ImagingSpatialObject,
     return sc_obj
 end
 
-function find_markers(sc_obj::Union{scRNAObject, VisiumObject, ImagingSpatialObject, CartanaObject, XeniumObject, scATACObject, MerfishObject, SlideseqObject, STARmapObject, seqFishObject}; cluster_1::Union{String, Nothing}=nothing, cluster_2::Union{String, Nothing}=nothing, 
+function find_markers(sc_obj::Union{scRNAObject, VisiumObject, ImagingSpatialObject, CartanaObject, XeniumObject, scATACObject, MerfishObject, SlideseqObject, STARmapObject, seqFishObject, StereoSeqObject}; cluster_1::Union{String, Nothing}=nothing, cluster_2::Union{String, Nothing}=nothing, 
     anno::Union{String, Symbol}="cluster", expr_cutoff=0.0, min_pct=0.1, p_cutoff = 0.05, only_pos = true)
     if isa(cluster_1, Nothing)
         error("Please provide the name of the cell cluster for which you wish to obtain the differential genes. The \"cluster_1\" parameter cannot be left blank.")
@@ -275,7 +275,7 @@ function find_markers(sc_obj::Union{scRNAObject, VisiumObject, ImagingSpatialObj
     return test_result
 end
 
-function find_all_markers(sc_obj::Union{scRNAObject, VisiumObject, ImagingSpatialObject, CartanaObject, XeniumObject, scATACObject,MerfishObject, SlideseqObject, STARmapObject, seqFishObject}; anno::Union{String, Symbol}="cluster", expr_cutoff=0.0, min_pct=0.1, p_cutoff = 0.05, only_pos = true)
+function find_all_markers(sc_obj::Union{scRNAObject, VisiumObject, ImagingSpatialObject, CartanaObject, XeniumObject, scATACObject,MerfishObject, SlideseqObject, STARmapObject, seqFishObject, StereoSeqObject}; anno::Union{String, Symbol}="cluster", expr_cutoff=0.0, min_pct=0.1, p_cutoff = 0.05, only_pos = true)
     if isa(anno, String)
         anno = Symbol(anno)
     end
diff --git a/src/scrna/utils.jl b/src/scrna/utils.jl
index df8485e..c216a67 100644
--- a/src/scrna/utils.jl
+++ b/src/scrna/utils.jl
@@ -1,8 +1,8 @@
 colSum(mtx::AbstractMatrix{<:Real}) = sum(mtx, dims=1)
 rowSum(mtx::AbstractMatrix{<:Real}) = sum(mtx, dims=2)
-rownames(sc_obj::Union{scRNAObject, VisiumObject, ImagingSpatialObject,  CartanaObject, XeniumObject, MerfishObject, SlideseqObject, seqFishObject, STARmapObject}) = sc_obj.rawCount.gene_name
+rownames(sc_obj::Union{scRNAObject, VisiumObject, ImagingSpatialObject,  CartanaObject, XeniumObject, MerfishObject, SlideseqObject, seqFishObject, STARmapObject, StereoSeqObject}) = sc_obj.rawCount.gene_name
 rownames(sc_obj::scATACObject) = sc_obj.rawCount.peak_name
-colnames(sc_obj::Union{scRNAObject, VisiumObject, ImagingSpatialObject, CartanaObject, XeniumObject, scATACObject, MerfishObject, SlideseqObject, seqFishObject, STARmapObject}) = sc_obj.rawCount.cell_name
+colnames(sc_obj::Union{scRNAObject, VisiumObject, ImagingSpatialObject, CartanaObject, XeniumObject, scATACObject, MerfishObject, SlideseqObject, seqFishObject, STARmapObject, StereoSeqObject}) = sc_obj.rawCount.cell_name
 rownames(ct_mat::AbstractCount) = ct_mat.gene_name
 colnames(ct_mat::AbstractCount) = ct_mat.cell_name
 
@@ -86,7 +86,7 @@ function subset_count(ct_obj::T;
     return new_obj
 end
 
-function extract_cluster_count(sc_obj::Union{scRNAObject, VisiumObject, ImagingSpatialObject, CartanaObject, XeniumObject, MerfishObject, SlideseqObject, seqFishObject, STARmapObject}, cl; count_type = "norm", anno = Union{String, Symbol}="cluster")
+function extract_cluster_count(sc_obj::Union{scRNAObject, VisiumObject, ImagingSpatialObject, CartanaObject, XeniumObject, MerfishObject, SlideseqObject, seqFishObject, STARmapObject, StereoSeqObject}, cl; count_type = "norm", anno = Union{String, Symbol}="cluster")
     df = sc_obj.clustData.clustering
     if isa(anno, String)
         anno = Symbol(anno)
@@ -150,12 +150,12 @@ function jitter(x)
     end
 end
 
-function variable_genes(sc_obj::Union{scRNAObject, VisiumObject, ImagingSpatialObject, CartanaObject, XeniumObject, MerfishObject, SlideseqObject, seqFishObject, STARmapObject})
+function variable_genes(sc_obj::Union{scRNAObject, VisiumObject, ImagingSpatialObject, CartanaObject, XeniumObject, MerfishObject, SlideseqObject, seqFishObject, STARmapObject, StereoSeqObject})
     vargenes = pbmc.varGene.var_gene
     return vargenes
 end
 
-function update_object(sp_obj::Union{scRNAObject, VisiumObject, ImagingSpatialObject, CartanaObject, XeniumObject, MerfishObject, SlideseqObject, seqFishObject, STARmapObject})
+function update_object(sp_obj::Union{scRNAObject, VisiumObject, ImagingSpatialObject, CartanaObject, XeniumObject, MerfishObject, SlideseqObject, seqFishObject, STARmapObject, StereoSeqObject})
     cells = colnames(sp_obj)
     genes = rownames(sp_obj)
     all_cells = sp_obj.metaData.Cell_id
@@ -189,7 +189,7 @@ function update_object(sp_obj::Union{scRNAObject, VisiumObject, ImagingSpatialOb
             sp_obj.clustData.adj_mat = sp_obj.clustData.adj_mat[check_cell, check_cell]
         end
     end
-    if isa(sp_obj, Union{ImagingSpatialObject, CartanaObject, XeniumObject, MerfishObject, seqFishObject, STARmapObject})
+    if isa(sp_obj, Union{ImagingSpatialObject, CartanaObject, XeniumObject, MerfishObject, seqFishObject, STARmapObject, StereoSeqObject})
         println("Updating spatial data...")
         sp_obj.spmetaData.cell = filter(:cell => ∈(cell_set), sp_obj.spmetaData.cell)
         sp_obj.spmetaData.molecule = filter(:cell => ∈(cell_set), sp_obj.spmetaData.molecule)
@@ -250,13 +250,13 @@ function update_object(sp_obj::Union{scRNAObject, VisiumObject, ImagingSpatialOb
     return sp_obj
 end
 
-function subset_object(sp_obj::Union{scRNAObject, VisiumObject, ImagingSpatialObject, CartanaObject, XeniumObject, MerfishObject, SlideseqObject, seqFishObject, STARmapObject}; cells = nothing, genes = nothing)
+function subset_object(sp_obj::Union{scRNAObject, VisiumObject, ImagingSpatialObject, CartanaObject, XeniumObject, MerfishObject, SlideseqObject, seqFishObject, STARmapObject, StereoSeqObject}; cells = nothing, genes = nothing)
     sp_obj.rawCount = subset_count(sp_obj.rawCount; genes = genes, cells = cells)
     sp_obj = update_object(sp_obj)
     return sp_obj
 end
 
-function check_dim(sp_obj::Union{scRNAObject, VisiumObject,ImagingSpatialObject, CartanaObject, XeniumObject, MerfishObject, seqFishObject, STARmapObject}, field_name::Union{Symbol, String})
+function check_dim(sp_obj::Union{scRNAObject, VisiumObject,ImagingSpatialObject, CartanaObject, XeniumObject, MerfishObject, seqFishObject, STARmapObject, StereoSeqObject}, field_name::Union{Symbol, String})
     if isa(field_name, String)
        field_name = Symbol(field_name)
     end
@@ -265,7 +265,7 @@ function check_dim(sp_obj::Union{scRNAObject, VisiumObject,ImagingSpatialObject,
     return check_length
    end
    
-function update_count(sp_obj::Union{scRNAObject, VisiumObject, ImagingSpatialObject, CartanaObject, XeniumObject, MerfishObject, seqFishObject, STARmapObject}, ct_name::Union{Symbol, String})
+function update_count(sp_obj::Union{scRNAObject, VisiumObject, ImagingSpatialObject, CartanaObject, XeniumObject, MerfishObject, seqFishObject, STARmapObject, StereoSeqObject}, ct_name::Union{Symbol, String})
     cell_id = colnames(sp_obj)
     gene_id = rownames(sp_obj)
     if isa(ct_name, String)
@@ -337,7 +337,7 @@ function sparse_r_to_jl(counts::RObject{S4Sxp})
     return julia_sparse_matrix
 end
 
-function annotate_cells(sp::Union{scRNAObject, VisiumObject, ImagingSpatialObject, CartanaObject, XeniumObject, MerfishObject, SlideseqObject, seqFishObject, STARmapObject}, cell_map::Dict; old_id_name::Union{String, Symbol}="cluster", new_id_name::Union{String, Symbol}="celltype")
+function annotate_cells(sp::Union{scRNAObject, VisiumObject, ImagingSpatialObject, CartanaObject, XeniumObject, MerfishObject, SlideseqObject, seqFishObject, STARmapObject, StereoSeqObject}, cell_map::Dict; old_id_name::Union{String, Symbol}="cluster", new_id_name::Union{String, Symbol}="celltype")
     sp.metaData = map_values(sp.metaData, old_id_name , new_id_name, collect(keys(cell_map)),  collect(values(cell_map)))
     if isa(sp, Union{CartanaObject, MerfishObject, XeniumObject, seqFishObject, STARmapObject})
         sp.spmetaData.cell[!, old_id_name] = sp.metaData.cell[!, old_id_name]
diff --git a/src/spatial/baysor_boundary.jl b/src/spatial/baysor_boundary.jl
new file mode 100644
index 0000000..d48ab11
--- /dev/null
+++ b/src/spatial/baysor_boundary.jl
@@ -0,0 +1,288 @@
+## All the functions below were from Baysor. Please visit the original repo to learn more: https://github.com/kharchenkolab/Baysor/tree/master
+
+function split(vector::T where T <: AbstractVector; n_parts::Int)
+    offset = ceil(Int, length(vector) / n_parts)
+    return [vector[n:min(n + offset - 1, length(vector))] for n in 1:offset:length(vector)]
+end
+
+function split(array::T where T <: AbstractVector{TV}, factor::T2 where T2 <: AbstractVector{<:Integer}; max_factor::Union{Int, Nothing}=nothing, drop_zero::Bool=false)::Array{Vector{TV}, 1} where TV
+    @assert length(array) == length(factor)
+    if max_factor === nothing
+        max_factor = maximum(factor)
+    end
+
+    splitted = [TV[] for i in 1:max_factor]
+    for i in 1:length(array)
+        if drop_zero && factor[i] == 0
+            continue
+        end
+
+        push!(splitted[factor[i]], array[i])
+    end
+
+    return splitted
+end
+
+split(array::UnitRange{Int64}, factor::Array{Int64,1}; kwargs...) = split(collect(array), factor; kwargs...)
+
+split_ids(factor::Array{Int, 1}; kwargs...) = split(1:length(factor), factor; kwargs...)
+
+split(df::DataFrame, factor::Symbol; kwargs...) = split(df, Array(df[!, factor]); kwargs...)
+split(df::DataFrame, factor::Vector{Int}; kwargs...) = [df[ids, :] for ids in split(1:size(df, 1), factor; kwargs...)]
+
+function grid_borders_per_label(grid_labels::Matrix{<:Unsigned})
+    d_cols = [-1 0 1 0]
+    d_rows = [0 1 0 -1]
+    borders_per_label = [Array{UInt32, 1}[] for i in 1:maximum(grid_labels)]
+
+    for row in 1:size(grid_labels, 1)
+        for col in 1:size(grid_labels, 2)
+            cur_label = grid_labels[row, col];
+            if cur_label == 0
+                continue
+            end
+
+            for (d_row, d_col) in zip(d_rows, d_cols)
+                n_row = row + d_row
+                n_col = col + d_col
+
+                if (n_row < 1) || (n_col < 1) || (n_row > size(grid_labels, 1)) || (n_col > size(grid_labels, 2)) ||
+                        (grid_labels[n_row, n_col] != cur_label)
+                    push!(borders_per_label[cur_label], [col, row])
+                    break
+                end
+            end
+        end
+    end
+
+    return borders_per_label
+end
+
+function border_mst(border_points::Array{<:Vector{<:Real},1}; min_edge_length::Float64=0.1, max_edge_length::Float64=1.5)
+    leafsize = min(size(border_points, 1), 10) # Some performance feature?
+    tree = KDTree(Matrix{Float64}(hcat(border_points...)), leafsize=leafsize);
+    src_verts, dst_verts, weights = Int[], Int[], Float64[]
+    edges = Set{Tuple{Int, Int}}()
+
+    n_neighbs = min(length(border_points), 5)
+    for (i, p) in enumerate(border_points)
+        inds, dists = knn(tree, p, n_neighbs, true) # 5 to get border for sure
+        filt_mask = (dists .> min_edge_length) .& (dists .< max_edge_length)
+
+        for (ind, dist) in zip(inds[filt_mask], dists[filt_mask])
+            edge = sort([i, ind])
+            edge = (edge...,)
+            if !(edge in edges)
+                push!(src_verts, edge[1])
+                push!(dst_verts, edge[2])
+                push!(weights, dist)
+                push!(edges, edge)
+            end
+        end
+    end
+
+    adj_mtx = sparse(src_verts, dst_verts, weights, size(border_points, 1), size(border_points, 1));
+    return Graphs.kruskal_mst(SimpleWeightedGraph(adj_mtx + adj_mtx'));
+end
+
+function find_longest_paths(edges::Array{<:Graphs.AbstractEdge, 1})::Array{Vector{Int}, 1}
+    if length(edges) == 0
+        return []
+    end
+
+    vert_counts = StatsBase.countmap(vcat(map(x -> [src(x), dst(x)], edges)...));
+    n_verts = maximum(keys(vert_counts));
+    leafs = collect(keys(vert_counts))[collect(values(vert_counts)) .== 1];
+
+    adj_mtx = sparse(src.(edges), dst.(edges), weight.(edges), n_verts, n_verts)
+    g_inc = SimpleWeightedGraph(adj_mtx + adj_mtx')
+
+    conn_comps = Graphs.connected_components(g_inc);
+    longest_paths = Array{Int, 1}[]
+    for vert_inds in conn_comps
+        if size(vert_inds, 1) < 3
+            continue
+        end
+
+        comp_leafs = intersect(vert_inds, leafs)
+        max_dist, max_source, max_target = -1, -1, -1
+        paths = []
+        for src_id in comp_leafs
+            paths = Graphs.dijkstra_shortest_paths(g_inc, src_id);
+            cur_target = comp_leafs[findmax(paths.dists[comp_leafs])[2]];
+            if isinf(paths.dists[cur_target])
+                @warn "Infinite distance for verteces $src_id and $cur_target"
+                continue
+            end
+
+            if paths.dists[cur_target] > max_dist
+                max_dist = paths.dists[cur_target]
+                max_target = cur_target
+                max_source = src_id
+            end
+        end
+
+        if max_source < 0 | max_target < 0
+            @warn "Source/target vertex id is 0"
+            continue
+        end
+
+        cur_target = max_target
+        path = Int[cur_target]
+        paths = Graphs.dijkstra_shortest_paths(g_inc, max_source);
+        if isinf(paths.dists[max_target])
+            @warn "Infinite distance for verteces $src_id and $cur_target"
+            continue
+        end
+        while cur_target != max_source
+            cur_target = paths.parents[cur_target]
+            push!(path, cur_target)
+            if size(path, 1) > 2 * paths.dists[max_target]
+                @warn "Can't build path"
+                break
+            end
+        end
+        push!(longest_paths, path)
+    end
+
+    return longest_paths
+end
+
+function order_points_to_polygon(vert_inds::Vector{Int}, border_coords::Matrix{T} where T <: Real; max_dev::T2 where T2 <: Real=10.0)
+    polygon_inds = Vector{Int}[]
+    while !isempty(vert_inds)
+        ci = vert_inds[1]
+        cur_ids = [ci]
+        kdt = KDTree(Float64.(border_coords[:, vert_inds]));
+        while true
+            nn_ids = setdiff(vert_inds[inrange(kdt, Float64.(border_coords[:, ci]), max_dev)], cur_ids)
+            if length(nn_ids) == 0
+                break
+            end
+
+            nn_ids = nn_ids[sortperm(vec(sum((border_coords[:, ci] .- border_coords[:, nn_ids]) .^ 2, dims=1)))]
+
+            ci = nn_ids[1]
+            push!(cur_ids, ci)
+        end
+
+        push!(polygon_inds, cur_ids)
+        vert_inds = setdiff(vert_inds, cur_ids)
+    end
+
+    return polygon_inds
+end
+
+function find_grid_point_labels_kde(pos_data::Matrix{Float64}, cell_labels::Vector{<:Integer}, min_x::Union{Vector{Float64}, Tuple{Float64, Float64}},
+        max_x::Union{Vector{Float64}, Tuple{Float64, Float64}}; grid_step::Float64, bandwidth::Float64, dens_threshold::Float64=1e-5, min_molecules_per_cell::Int=3,
+        verbose::Bool=false)::Matrix{<:Unsigned}  where T <: Real
+    coords_per_label = [pos_data[:, ids] for ids in split_ids(cell_labels .+ 1)];
+    filt_mask = (size.(coords_per_label, 2) .>= min_molecules_per_cell)
+    coords_per_label = coords_per_label[filt_mask];
+    cell_labels = (0:maximum(cell_labels))[filt_mask]
+
+    if length(coords_per_label) == 0
+        return Matrix{Int}(undef,0,0)
+    end
+
+    FFTW.set_num_threads(1); # see https://github.com/JuliaStats/KernelDensity.jl/issues/80
+
+    xw, yw = ceil.(Int, (max_x .- min_x) ./ grid_step)
+    label_mat = nothing
+    if maximum(cell_labels) < (2 ^ 16) - 1
+        label_mat = zeros(UInt16, yw, xw)
+    else
+        label_mat = zeros(UInt32, yw, xw)
+    end
+    dens_mat = zeros(Float16, yw, xw)
+
+    p = verbose ? Progress(length(coords_per_label)) : nothing
+    for (ci, lab) in zip(1:length(coords_per_label), cell_labels)
+        c_bandwidth = (ci == 1) ? bandwidth / 1.5 : bandwidth
+        coords = coords_per_label[ci]
+        sxi = max(floor(Int, (minimum(coords[1,:]) - min_x[1] - 3 * c_bandwidth) / grid_step), 1);
+        ex = maximum(coords[1,:]) + 3 * c_bandwidth;
+
+        syi = max(floor(Int, (minimum(coords[2,:]) - min_x[2] - 3 * c_bandwidth) / grid_step), 1);
+        ey = maximum(coords[2,:]) + 3 * c_bandwidth;
+
+        cxs = (sxi * grid_step + min_x[1]):grid_step:ex;
+        cys = (syi * grid_step + min_x[2]):grid_step:ey;
+
+        dens = kde((coords[2,:], coords[1,:]), (cys, cxs), bandwidth=(c_bandwidth, c_bandwidth)).density
+
+        for dyi in 1:length(cys)
+            yi = dyi + syi
+            if yi > size(dens_mat, 1)
+                break
+            end
+
+            for dxi in 1:length(cxs)
+                xi = dxi + sxi
+                if xi > size(dens_mat, 2)
+                    break
+                end
+                d = dens[dyi, dxi] * size(coords, 2)
+                if (d > dens_threshold) && (d > dens_mat[yi, xi])
+                    dens_mat[yi, xi] = d
+                    label_mat[yi, xi] = lab
+                end
+            end
+        end
+
+        if verbose
+            next!(p)
+        end
+    end
+
+    return ImageMorphology.closing(label_mat)
+end
+
+function extract_polygons_from_label_grid(grid_labels::Matrix{<:Unsigned}; min_border_length::Int=3, shape_method::Symbol=:path, max_dev::Float64=10.0,
+        exclude_labels::Vector{Int}=Int[], offset::Vector{Float64}=[0.0, 0.0], grid_step::Float64=1.0)::Array{Matrix{Float64}, 1}
+    borders_per_label = grid_borders_per_label(grid_labels);
+    if !isempty(exclude_labels)
+        borders_per_label = borders_per_label[setdiff(1:length(borders_per_label), exclude_labels)]
+    end
+
+    borders_per_label = borders_per_label[length.(borders_per_label) .>= min_border_length]
+
+    paths = nothing
+    edges_per_label = border_mst.(borders_per_label)
+    if shape_method == :path
+        paths = find_longest_paths.(edges_per_label);
+    elseif shape_method == :order
+        conn_comps = [Graphs.connected_components(Graphs.SimpleGraphFromIterator(Graphs.SimpleGraphs.SimpleEdge.(src.(edges), dst.(edges)))) for edges in edges_per_label]
+        paths = [vcat(order_points_to_polygon.(cc, Ref(hcat(borders...)), max_dev=max_dev)...) for (borders, cc) in zip(borders_per_label, conn_comps)];
+    else
+        error("Unknown shape_method: $shape_method")
+    end
+
+    paths = [ps[length.(ps) .> min_border_length] for ps in paths]
+    polygons = vcat([[borders[p] for p in cur_paths] for (borders, cur_paths) in zip(borders_per_label, paths)]...);
+    polygons = [copy((hcat(coords...) .* grid_step .+ offset)') for coords in polygons]
+    return [vcat(cp, cp[1,:]') for cp in polygons]
+end
+
+function baysor_boundary_polygons(molecules::DataFrame, cell_labels::Vector{<:Integer}; min_x::Union{Array{Float64}, Nothing}=nothing, max_x::Union{Array{Float64}, Nothing}=nothing,
+                           grid_step::Float64=5.0, min_border_length::Int=3, shape_method::Symbol=:path, max_dev::Float64=10.0, x_col=:x, y_col=:y,
+                           bandwidth::Float64=(grid_step / 2), exclude_labels::Vector{Int}=Int[], kwargs...)::Array{Matrix{Float64}, 1}
+    pos_data = copy(Matrix{Float64}(molecules[:, [x_col, y_col]])')
+    if min_x === nothing
+        min_x = vec(mapslices(minimum, pos_data, dims=2))
+    end
+
+    if max_x === nothing
+        max_x = vec(mapslices(maximum, pos_data, dims=2))
+    end
+
+    grid_labels_plane = find_grid_point_labels_kde(pos_data, cell_labels, min_x, max_x; grid_step=grid_step, bandwidth=bandwidth, kwargs...)
+    GC.gc()
+
+    if length(grid_labels_plane) == 0
+        return Matrix{Float64}[]
+    end
+
+    return extract_polygons_from_label_grid(grid_labels_plane; min_border_length=min_border_length, shape_method=shape_method, max_dev=max_dev,
+        exclude_labels=exclude_labels, offset=min_x, grid_step=grid_step)
+end
diff --git a/src/spatial/sp_analysis.jl b/src/spatial/sp_analysis.jl
index 4fd4783..1df9e30 100644
--- a/src/spatial/sp_analysis.jl
+++ b/src/spatial/sp_analysis.jl
@@ -51,7 +51,7 @@ function run_SpaGCN(sp::AbstractSpaObj, count_path::String, python_path::String;
             return sp
 end
 
-function run_tangram(sp_obj::Union{ImagingSpatialObject, CartanaObject, XeniumObject,MerfishObject, STARmapObject, seqFishObject}, 
+function run_tangram(sp_obj::Union{ImagingSpatialObject, CartanaObject, XeniumObject,MerfishObject, STARmapObject, seqFishObject, StereoSeqObject}, 
     sc_obj::scRNAObject; gene_list::Union{String, Vector{String}, Nothing}=nothing,
     density_prior="uniform",num_epochs=100,device="cpu")
     sc = pyimport("scanpy")
@@ -106,7 +106,7 @@ function run_tangram(sp_obj::Union{ImagingSpatialObject, CartanaObject, XeniumOb
     return sp_obj
 end
 
-function run_spaGE(sp::Union{ImagingSpatialObject, CartanaObject, XeniumObject,MerfishObject, STARmapObject, seqFishObject}, 
+function run_spaGE(sp::Union{ImagingSpatialObject, CartanaObject, XeniumObject,MerfishObject, STARmapObject, seqFishObject, StereoSeqObject}, 
     sc::scRNAObject, spaGE_path::String; gene_list::Union{String, Vector{String}, Nothing}=nothing, npv::Int64=30)
     if isa(gene_list, String)
         gene_list = [gene_list]
@@ -155,7 +155,7 @@ function run_spaGE(sp::Union{ImagingSpatialObject, CartanaObject, XeniumObject,M
     return sp
 end
 
-function run_gimVI(sp_obj::Union{ImagingSpatialObject, CartanaObject, XeniumObject,MerfishObject, STARmapObject, seqFishObject}, 
+function run_gimVI(sp_obj::Union{ImagingSpatialObject, CartanaObject, XeniumObject,MerfishObject, STARmapObject, seqFishObject, StereoSeqObject}, 
     sc_obj::scRNAObject; gene_list::Union{String, Vector{String}, Nothing}=nothing, epochs::Int64=200)
     if isa(gene_list, String)
         gene_list = [gene_list]
diff --git a/src/spatial/sp_objects.jl b/src/spatial/sp_objects.jl
index 58a6d57..e3bc6aa 100644
--- a/src/spatial/sp_objects.jl
+++ b/src/spatial/sp_objects.jl
@@ -542,3 +542,71 @@ mutable struct seqFishObject <: AbstractImagingObj
         println("seqFishObject was successfully created!")
     end
 end
+
+mutable struct SpaImageObj <: AbstractImagingObj
+    histoImage::Union{Matrix{RGB{N0f8}},Matrix{Gray{N0f8}}, Nothing}
+    ifImage::Union{Matrix{RGB{N0f8}},Matrix{Gray{N0f8}}, Nothing}
+    dapiImage::Union{Matrix{RGB{N0f8}},Matrix{Gray{N0f8}}, Nothing}
+    SpaImageObj(histoImage, ifImage, dapiImage) = new(histoImage, ifImage, dapiImage)
+end
+
+mutable struct StereoSeqObject <: AbstractImagingObj
+    rawCount::Union{RawCountObject, Nothing}
+    normCount::Union{NormCountObject, Nothing}
+    scaleCount::Union{ScaleCountObject, Nothing}
+    metaData::Union{DataFrame, Nothing}
+    spmetaData::Union{SpaMetaObj, Nothing}
+    varGene::Union{VariableGeneObject, Nothing}
+    dimReduction::Union{ReductionObject, Nothing}
+    clustData::Union{ClusteringObject, Nothing}
+    polyCount::Union{RawCountObject, Nothing}
+    polynormCount::Union{NormCountObject, Nothing}
+    coordData::Union{SpaCoordObj, Nothing}
+    imputeData::Union{SpaImputeObj, Nothing}
+    imgObject::Union{SpaImageObj, Nothing}
+    polygonData::Union{Vector{Matrix{Float64}}, Nothing}
+    function StereoSeqObject(molecule_data::DataFrame, cell_data::DataFrame, counts::RawCountObject; 
+        prefix::Union{String, Nothing}=nothing, postfix::Union{String, Nothing}=nothing, meta_data::Union{DataFrame, Nothing} = nothing,
+        min_gene::Int64=0, min_cell::Int64=0, x_col::Union{String, Symbol} = "x", 
+        y_col::Union{String, Symbol} = "y", cell_col::Union{String, Symbol} = "cell")
+        if isa(prefix, String)
+            println("Adding prefix " * prefix * " to all cells...")
+            counts.cell_name = prefix * "_" .* counts.cell_name
+            molecule_data[!, cell_col] = prefix * "_" .* molecule_data[!, cell_col]
+            cell_data[!, cell_col] = prefix * "_" .* cell_data[!, cell_col]
+        end
+        if isa(postfix, String)
+            println("Adding postfix " * postfix * " to all cells...")
+            counts.cell_name = counts.cell_name .* "_" .* postfix
+            molecule_data[!, cell_col] = molecule_data[!, cell_col] .* "_" .* postfix
+            cell_data[!, cell_col] = cell_data[!, cell_col] .* "_" .* postfix
+        end
+        count_mat = counts.count_mtx
+        gene_name = counts.gene_name
+        cell_name = counts.cell_name
+        if min_gene > 0 || min_cell > 0
+            count_mat, gene_name, cell_name = subset_matrix(count_mat, gene_name, cell_name, min_gene, min_cell)
+            cell_check = check_vec(cell_name, cell_data[!, cell_col])
+            cell_data = cell_data[cell_check, :]
+            mol_check = check_vec(cell_name, molecule_data[!, cell_col])
+            molecule_data = molecule_data[mol_check, :]
+        end
+        if isa(meta_data, Nothing)
+            nFeatures = vec(colSum(count_mat))
+            nGenes = vec(sum(x->x>0, count_mat, dims=1))
+            meta_data = DataFrame(Cell_id = cell_name, nFeatures=nFeatures, nGenes = nGenes)
+        end
+        counts = RawCountObject(count_mat, cell_name, gene_name)
+        spObj=new(counts)
+        meta = SpaMetaObj(cell_data, molecule_data, nothing)
+        spObj.spmetaData = meta
+        cell_coord = cell_data[!, [x_col, y_col]]
+        mol_coord = molecule_data[!, [x_col, y_col]]
+        coord = SpaCoordObj(cell_coord, mol_coord, nothing, nothing)
+        spObj.coordData = coord
+        spObj.metaData = meta_data
+        spObj.polygonData = nothing
+        return spObj
+        println("StereoSeqObject was successfully created!")
+    end
+end
\ No newline at end of file
diff --git a/src/spatial/sp_plots.jl b/src/spatial/sp_plots.jl
index 3887407..876d047 100644
--- a/src/spatial/sp_plots.jl
+++ b/src/spatial/sp_plots.jl
@@ -1,5 +1,5 @@
 
-function sp_dot_plot(sp::Union{ImagingSpatialObject, CartanaObject, VisiumObject, XeniumObject, MerfishObject, SlideseqObject, STARmapObject, seqFishObject}, genes::Union{Vector, String},
+function sp_dot_plot(sp::Union{ImagingSpatialObject, CartanaObject, VisiumObject, XeniumObject, MerfishObject, SlideseqObject, STARmapObject, seqFishObject, StereoSeqObject}, genes::Union{Vector, String},
     cluster::Union{Symbol, String};expr_cutoff::Union{Float64, Int64}=0, split_by::Union{String, Nothing}=nothing,
     x_title="Gene",y_title="Cell type", use_imputed= false, imp_type = "SpaGE", cell_order::Union{Vector, String, Nothing}=nothing,
     fontsize::Int64=12, color_scheme::String="yelloworangered",reverse_color::Bool=false,
@@ -87,7 +87,7 @@ function sp_dot_plot(sp::Union{ImagingSpatialObject, CartanaObject, VisiumObject
     return p
 end
 
-function sp_feature_plot(sp::Union{ImagingSpatialObject, CartanaObject, VisiumObject, XeniumObject,MerfishObject, SlideseqObject, STARmapObject, seqFishObject}, gene_list::Union{String, Vector{String}, Tuple{String}}; layer::String = "cells", x_col::Union{String, Symbol}="x",
+function sp_feature_plot(sp::Union{ImagingSpatialObject, CartanaObject, VisiumObject, XeniumObject,MerfishObject, SlideseqObject, STARmapObject, seqFishObject, StereoSeqObject}, gene_list::Union{String, Vector{String}, Tuple{String}}; layer::String = "cells", x_col::Union{String, Symbol}="x",
     y_col::Union{String, Symbol}="y", cell_col = "cell", x_lims=nothing, y_lims=nothing, marker_size=2, order::Bool=true, scale::Bool = false,titlesize::Int64=24, 
     height::Real = 500, width::Real = 500, combine = true, img_res::String = "low",  adjust_contrast::Real = 1.0, adjust_brightness::Real = 0.3, use_imputed=false, imp_type::Union{String, Nothing} = nothing,
     color_keys=["gray94","orange","red3"], gene_colors = nothing, alpha = [1.0,1.0], clip = 0, legend_fontsize = 10, do_legend=false, legend_size = 10, bg_color = "white")
@@ -226,8 +226,16 @@ function sp_feature_plot(sp::Union{ImagingSpatialObject, CartanaObject, VisiumOb
                 df_plt[!, x_col] = df_plt[!, x_col] .- x_lims[1]
                 df_plt[!, y_col] = df_plt[!, y_col] .- y_lims[1]
             else
-                MK.xlims!(MK.current_axis(), x_lims1)
-                MK.ylims!(MK.current_axis(), y_lims1)
+                if isa(x_lims, Nothing)
+                    MK.xlims!(MK.current_axis(), x_lims1)
+                else
+                    MK.xlims!(MK.current_axis(), x_lims)
+                end
+                if isa(y_lims, Nothing)
+                    MK.ylims!(MK.current_axis(), y_lims1)
+                else
+                    MK.ylims!(MK.current_axis(), y_lims)
+                end                    
             end
             MK.scatter!(ax1, df_plt[!, x_col], df_plt[!, y_col]; color = df_plt.plt_color, strokewidth = 0, markersize = marker_size)
             MK.Colorbar(fig[n_row,n_col2], label = "", colormap = c_map, width=10, limits = (0, maximum(gene_expr)))
@@ -334,7 +342,7 @@ function sp_feature_plot(sp::Union{ImagingSpatialObject, CartanaObject, VisiumOb
     end
 end
 
-function plot_gene_polygons(sp::Union{ImagingSpatialObject, CartanaObject,XeniumObject,MerfishObject, STARmapObject, seqFishObject}, gene_list::Union{String, Vector{String}, Tuple{String}};
+function plot_gene_polygons(sp::Union{ImagingSpatialObject, CartanaObject,XeniumObject,MerfishObject, STARmapObject, seqFishObject, StereoSeqObject}, gene_list::Union{String, Vector{String}, Tuple{String}};
     color_keys::Union{Vector{String}, Tuple{String,String,String}}=["gray94","lemonchiffon","orange","red3"],
     x_lims=nothing, y_lims=nothing,width=900,height=1000,stroke_width=0,stroke_color="black",
     titlesize::Int64=24, scale::Bool = false, bg_color = "white"
@@ -393,7 +401,7 @@ function plot_gene_polygons(sp::Union{ImagingSpatialObject, CartanaObject,Xenium
         MK.current_figure()
 end
 
-function plot_cell_polygons(sp::Union{ImagingSpatialObject, CartanaObject,XeniumObject,MerfishObject, STARmapObject, seqFishObject}, anno;
+function plot_cell_polygons(sp::Union{ImagingSpatialObject, CartanaObject,XeniumObject,MerfishObject, STARmapObject, seqFishObject, StereoSeqObject}, anno;
     cell_types::Union{String, Int64, Vector, Tuple, Nothing}=nothing, cell_colors::Union{Nothing, String, Vector, Tuple} = nothing,
     pt_bg_color = "gray90", x_lims=nothing, y_lims=nothing,width = 500, height = 500,stroke_width=0, stroke_color="black",
     legend_fontsize = 10, do_legend=false, legend_size = 10 , bg_color = "white"
@@ -459,13 +467,15 @@ function plot_cell_polygons(sp::Union{ImagingSpatialObject, CartanaObject,Xenium
     return MK.current_figure()
 end
 
-function plot_transcript_polygons(sp::Union{ImagingSpatialObject, CartanaObject,XeniumObject,MerfishObject, STARmapObject, seqFishObject}; 
+function plot_transcript_polygons(sp::Union{ImagingSpatialObject, CartanaObject,XeniumObject,MerfishObject, STARmapObject, seqFishObject, StereoSeqObject}; 
     gene_list::Union{Vector, Symbol, String, Nothing}=nothing, 
     gene_colors::Union{Vector, Symbol, String, Nothing}=nothing, 
     pt_bg_color::Union{Vector, Symbol, String}="gray95",
     stroke_color::Union{Vector, Symbol, String}="black",
     marker_size = 2, stroke_width=0.1, bg_color="white",
-    canvas_size=(900,1000),x_lims=nothing, y_lims=nothing, 
+    width = 900,
+    height = 1000,
+    x_lims=nothing, y_lims=nothing, 
     anno::Union{Symbol, String, Nothing}=nothing,
     ann_colors::Union{Nothing, Dict}=nothing
 )
@@ -496,10 +506,16 @@ function plot_transcript_polygons(sp::Union{ImagingSpatialObject, CartanaObject,
     if isa(anno, String)
         anno = Symbol(anno)
     end
-    anno_df=sp.spmetaData.polygon
+    select_fov = filter([:x, :y] => (x, y) -> x_lims[1] < x < x_lims[2] && y_lims[1] < y < y_lims[2], sp.spmetaData.cell)
+    cell_set = Set(select_fov.cell)
+    anno_df = filter(:mapped_cell => ∈(cell_set), sp.spmetaData.polygon)
+    polygon_num = anno_df.polygon_number
+    polygons = sp.polygonData
+    polygons = polygons[polygon_num]
+    df_spatial = filter(:cell => ∈(cell_set), df_spatial)
     if isa(ann_colors, Nothing)
         cell_anno=unique(anno_df[!,anno])
-        c_map=Colors.distinguishable_colors(length(cell_anno), Colors.colorant"#007a10", lchoices=range(20, stop=70, length=15))
+        c_map=Colors.distinguishable_colors(length(cell_anno), parse(Colors.Colorant, pt_bg_color), lchoices=range(20, stop=70, length=15))
         c_map = "#" .* hex.(c_map)
         ann_colors=Dict(cell_anno .=> c_map)
     end
@@ -509,7 +525,7 @@ function plot_transcript_polygons(sp::Union{ImagingSpatialObject, CartanaObject,
     colors1 = df_spatial1[!,:new_color]
     df_spatial2=filter(:gene => ∈(Set(gene_list)), df_spatial)
     colors2 = df_spatial2[!,:new_color]
-    fig = MK.Figure(resolution=canvas_size)
+    fig = MK.Figure(resolution=(width, height))
     fig[1, 1] = MK.Axis(fig; backgroundcolor = bg_color, xticklabelsize=12, yticklabelsize=12, 
         xticksvisible=false, xticklabelsvisible=false, yticksvisible=false, 
         yticklabelsvisible=false, xgridvisible = false, ygridvisible = false)
@@ -524,7 +540,7 @@ function plot_transcript_polygons(sp::Union{ImagingSpatialObject, CartanaObject,
     return MK.current_figure()
 end
 
-function sp_dim_plot(sp::Union{ImagingSpatialObject, CartanaObject, VisiumObject,XeniumObject,MerfishObject, SlideseqObject, STARmapObject, seqFishObject}, anno::Union{Symbol, String}; 
+function sp_dim_plot(sp::Union{ImagingSpatialObject, CartanaObject, VisiumObject,XeniumObject,MerfishObject, SlideseqObject, STARmapObject, seqFishObject, StereoSeqObject}, anno::Union{Symbol, String}; 
     anno_color::Union{Nothing, Dict} = nothing, x_col::String = "x", y_col::String = "y", cell_order::Union{Vector{String}, Nothing}=nothing,
     x_lims=nothing, y_lims=nothing,canvas_size=(900,1000),stroke_width=0.5,stroke_color=:transparent,  bg_color=:white,
         marker_size=2, label_size=50, label_color="black", label_offset=(0,0), do_label=false, do_legend=true, alpha::Real = 1,
@@ -641,7 +657,7 @@ function sp_dim_plot(sp::Union{ImagingSpatialObject, CartanaObject, VisiumObject
     return MK.current_figure()
 end
 
-function sp_highlight_cells(sp::Union{ImagingSpatialObject, CartanaObject, VisiumObject,XeniumObject,MerfishObject, SlideseqObject, STARmapObject, seqFishObject}, cell_hightlight::String, anno::Union{String,Symbol};
+function sp_highlight_cells(sp::Union{ImagingSpatialObject, CartanaObject, VisiumObject,XeniumObject,MerfishObject, SlideseqObject, STARmapObject, seqFishObject, StereoSeqObject}, cell_hightlight::String, anno::Union{String,Symbol};
     canvas_size=(900,1000),stroke_width::Float64=0.1, stroke_color="black", cell_color::String="red",
     marker_size=2,x_lims=nothing, y_lims=nothing, bg_color = "white")
     coord_cells=deepcopy(sp.spmetaData.cell)
@@ -665,7 +681,7 @@ function sp_highlight_cells(sp::Union{ImagingSpatialObject, CartanaObject, Visiu
     MK.current_figure()
 end
 
-function sp_gene_rank(sp::Union{ImagingSpatialObject, CartanaObject, XeniumObject, MerfishObject, STARmapObject, seqFishObject}, celltype::String, cluster::String; num_gene::Int64=20)
+function sp_gene_rank(sp::Union{ImagingSpatialObject, CartanaObject, XeniumObject, MerfishObject, STARmapObject, seqFishObject, StereoSeqObject}, celltype::String, cluster::String; num_gene::Int64=20)
     gene_list=unique(sp.spmetaData.molecule.gene)
     all_df=DataFrame()
     for (i, gene) in enumerate(gene_list)
@@ -690,7 +706,7 @@ function sp_gene_rank(sp::Union{ImagingSpatialObject, CartanaObject, XeniumObjec
         y={:rank,axis={title="Ranking", grid=false}})
 end
 
-function sp_feature_plot_group(sp_list::Union{ Vector{ImagingSpatialObject}, Vector{CartanaObject}, Vector{XeniumObject}, Vector{VisiumObject}, Vector{MerfishObject}, Vector{SlideseqObject}, Vector{seqFishObject}, Vector{STARmapObject}}, genes::Vector{String};
+function sp_feature_plot_group(sp_list::Union{ Vector{ImagingSpatialObject}, Vector{CartanaObject}, Vector{XeniumObject}, Vector{VisiumObject}, Vector{MerfishObject}, Vector{SlideseqObject}, Vector{seqFishObject}, Vector{STARmapObject}, Vector{StereoSeqObject}}, genes::Vector{String};
     x_col::Union{String, Symbol}="x",y_col::Union{String, Symbol}="y", alpha = [1.0,1.0], clip = 0,
     marker_size = 2, order=false, use_imputed = false,  imp_type::Union{String, Nothing}=nothing,
     height::Real = 500, width::Real = 500, titlesize::Int64 = 24, labels=nothing, img_res="low",
@@ -840,7 +856,7 @@ function sp_feature_plot_group(sp_list::Union{ Vector{ImagingSpatialObject}, Vec
     MK.current_figure()
 end
 
-function plot_fov(sp::Union{ImagingSpatialObject, CartanaObject,XeniumObject,MerfishObject, STARmapObject, seqFishObject}, n_fields_x::Int64, n_fields_y::Int64; 
+function plot_fov(sp::Union{ImagingSpatialObject, CartanaObject,XeniumObject,MerfishObject, STARmapObject, seqFishObject, StereoSeqObject}, n_fields_x::Int64, n_fields_y::Int64; 
     x_col::Union{String, Symbol}="x", y_col::Union{String, Symbol}="y", group_label::Union{Nothing, String}=nothing, 
     width=4000, height=4000, cell_highlight::Union{Nothing, String, Number}=nothing, shield::Bool= false, marker_size::Union{Int64, Float64}=10)
     df = sp.spmetaData.cell
@@ -886,7 +902,7 @@ function plot_fov(sp::Union{ImagingSpatialObject, CartanaObject,XeniumObject,Mer
     MK.current_figure()
 end
 
-function plot_point(sp::Union{ImagingSpatialObject, CartanaObject, VisiumObject, XeniumObject,MerfishObject, SlideseqObject, STARmapObject, seqFishObject}, pt::Vector{Float64}; 
+function plot_point(sp::Union{ImagingSpatialObject, CartanaObject, VisiumObject, XeniumObject,MerfishObject, SlideseqObject, STARmapObject, seqFishObject, StereoSeqObject}, pt::Vector{Float64}; 
     canvas_size=(4000,4000),marker_size=60, text_size=100, 
     pt_color="red", text_color="blue", label="point")
     df = sp.spmetaData.cell
@@ -905,7 +921,7 @@ function plot_point(sp::Union{ImagingSpatialObject, CartanaObject, VisiumObject,
     MK.current_figure()
 end
 
-function plot_depth(sp::Union{ImagingSpatialObject, CartanaObject, VisiumObject, XeniumObject,MerfishObject, SlideseqObject}; celltype::Union{String, Symbol} = :celltype,
+function plot_depth(sp::Union{ImagingSpatialObject, CartanaObject, VisiumObject, XeniumObject,MerfishObject, SlideseqObject, seqFishObject, StereoSeqObject}; celltype::Union{String, Symbol} = :celltype,
     cmap=nothing, cell_select=nothing, fontsize=16, scale=0.8, markers=nothing)
         cells=sp.spmetaData.cell
         celltypes=cell_select
@@ -958,7 +974,7 @@ function plot_depth(sp::Union{ImagingSpatialObject, CartanaObject, VisiumObject,
         MK.current_figure()
 end
 
-function plot_depth_animation(sp::Union{ImagingSpatialObject, CartanaObject, VisiumObject, XeniumObject,MerfishObject, SlideseqObject, STARmapObject, seqFishObject}, celltypes::Vector{String}, markers::Vector{String}; 
+function plot_depth_animation(sp::Union{ImagingSpatialObject, CartanaObject, VisiumObject, XeniumObject,MerfishObject, SlideseqObject, STARmapObject, seqFishObject, StereoSeqObject}, celltypes::Vector{String}, markers::Vector{String}; 
     group_label="celltype",gene_label="gene", bg_color="gray94",fontsize=16, scale=0.8, width=1800,height=600, file_name="animation.gif", framerate=30,
     titles = ["Cells colored by kidney depth","Cell distribution from cortex to papilla","Transcript distribution from cortex to papilla"])
     cells = sp.spmetaData.cell
@@ -1024,7 +1040,7 @@ function plot_depth_animation(sp::Union{ImagingSpatialObject, CartanaObject, Vis
         end
 end
 
-function plot_gene_depth(sp::Union{ImagingSpatialObject, CartanaObject, VisiumObject, XeniumObject,MerfishObject, SlideseqObject, STARmapObject, seqFishObject}, gene::String;
+function plot_gene_depth(sp::Union{ImagingSpatialObject, CartanaObject, VisiumObject, XeniumObject,MerfishObject, SlideseqObject, STARmapObject, seqFishObject, StereoSeqObject}, gene::String;
     c_map::Union{String, Symbol, Nothing}=nothing, cell_col="cell",
     canvas_size =(1200,300),marker_size=4,
     stroke_width=0.5,stroke_color="gray94",
@@ -1058,7 +1074,7 @@ function plot_gene_depth(sp::Union{ImagingSpatialObject, CartanaObject, VisiumOb
     MK.current_figure()
 end
 
-function PlotInteractive(sp::Union{CartanaObject, VisiumObject, XeniumObject,MerfishObject, SlideseqObject}; layer::String = "cells", marker_color::Union{Symbol, String}="black", marker_size=3, plot_mode="markers")
+function PlotInteractive(sp::Union{CartanaObject, VisiumObject, XeniumObject,MerfishObject, SlideseqObject, seqFishObject, StereoSeqObject}; layer::String = "cells", marker_color::Union{Symbol, String}="black", marker_size=3, plot_mode="markers")
     if layer === "cells"
         cells=sp.spmetaData.cell
         plyjs.plot(plyjs.scatter(x=cells.x, y=cells.y, mode=plot_mode, marker=plyjs.attr(size=marker_size, color=marker_color)))
@@ -1070,7 +1086,7 @@ function PlotInteractive(sp::Union{CartanaObject, VisiumObject, XeniumObject,Mer
     end
 end
 
-function plot_transcript_dapi(sp::Union{ImagingSpatialObject, CartanaObject,XeniumObject,MerfishObject, STARmapObject, seqFishObject}, fov::Int64, n_fields_x::Int64, 
+function plot_transcript_dapi(sp::Union{ImagingSpatialObject, CartanaObject,XeniumObject,MerfishObject, STARmapObject, seqFishObject, StereoSeqObject}, fov::Int64, n_fields_x::Int64, 
     n_fields_y::Int64; noise_ann=nothing,annotation=:cell,
     is_noise=nothing, draw_poly=false, marker_size=3)
     selected_view = subset_fov(sp, fov, n_fields_x, n_fields_y)
@@ -1112,7 +1128,7 @@ function plot_transcript_dapi(sp::Union{ImagingSpatialObject, CartanaObject,Xeni
     MK.current_figure()
 end
 
-function compare_gene_imputation(sp1::Union{ImagingSpatialObject, CartanaObject,XeniumObject, MerfishObject, STARmapObject, seqFishObject},sp2::Union{CartanaObject,XeniumObject,MerfishObject, STARmapObject, seqFishObject}, gene_list::Union{Vector, String},
+function compare_gene_imputation(sp1::Union{ImagingSpatialObject, CartanaObject,XeniumObject, MerfishObject, STARmapObject, seqFishObject, StereoSeqObject},sp2::Union{CartanaObject,XeniumObject,MerfishObject, STARmapObject, seqFishObject, StereoSeqObject}, gene_list::Union{Vector, String},
     cluster::Union{Symbol, String}; sp1_name::String ="sp1", sp2_name::String="sp2",
     assay_use::String="measured",expr_cutoff::Union{Float64, Int64}=0, legend_min::Union{Float64, Int64}=0, legend_max::Union{Float64, Int64}=1, 
     x_title="Gene",y_title="Cell type", cell_order::Union{Vector, String, Nothing}=nothing,
@@ -1179,7 +1195,7 @@ function compare_gene_imputation(sp1::Union{ImagingSpatialObject, CartanaObject,
         return p
 end
 
-function plot_heatmap(sp::Union{ImagingSpatialObject, CartanaObject,XeniumObject, MerfishObject, STARmapObject, seqFishObject}, gene_list::Union{Vector, String},
+function plot_heatmap(sp::Union{ImagingSpatialObject, CartanaObject,XeniumObject, MerfishObject, STARmapObject, seqFishObject, StereoSeqObject}, gene_list::Union{Vector, String},
     cluster::Union{Symbol, String};assay_use::String="measured",expr_cutoff::Union{Float64, Int64}=0, split_by::Union{String, Nothing}=nothing,
     x_title="Gene",y_title="Cell type", cell_order::Union{Vector, String, Nothing}=nothing,imp_type="SpaGE",
     fontsize::Int64=12, color_keys=["white", "ivory","gold","orange","tomato","red"],reverse_color::Bool=false,scale::Bool=false,
@@ -1264,7 +1280,7 @@ function plot_heatmap(sp::Union{ImagingSpatialObject, CartanaObject,XeniumObject
         return p
 end
 
-function overlay_visium_cartana(vs::Union{VisiumObject, SlideseqObject}, sp::Union{ImagingSpatialObject, CartanaObject,XeniumObject,MerfishObject, STARmapObject, seqFishObject}; vs_x = "new_x", vs_y = "new_y", 
+function overlay_visium_cartana(vs::Union{VisiumObject, SlideseqObject}, sp::Union{ImagingSpatialObject, CartanaObject,XeniumObject,MerfishObject, STARmapObject, seqFishObject, StereoSeqObject}; vs_x = "new_x", vs_y = "new_y", 
     sp_x = "new_x", sp_y = "new_y", vs_color=:red, sp_color=:blue, vs_markersize=7, 
     sp_markersize=2, vs_title="Visium", sp_title="Cartana")
     cartana_df = deepcopy(sp.spmetaData.cell)
@@ -1286,7 +1302,7 @@ function overlay_visium_cartana(vs::Union{VisiumObject, SlideseqObject}, sp::Uni
     MK.current_figure()
 end
 
-function overlay_visium_cartana_gene(vs::Union{VisiumObject, SlideseqObject}, sp::Union{ImagingSpatialObject, CartanaObject,XeniumObject, MerfishObject, STARmapObject, seqFishObject}, gene; vs_x="new_x", vs_y="new_y", sp_x="new_x", sp_y="new_y",
+function overlay_visium_cartana_gene(vs::Union{VisiumObject, SlideseqObject}, sp::Union{ImagingSpatialObject, CartanaObject,XeniumObject, MerfishObject, STARmapObject, seqFishObject, StereoSeqObject}, gene; vs_x="new_x", vs_y="new_y", sp_x="new_x", sp_y="new_y",
     vs_color=:red, sp_color=:blue, vs_markersize=7, canvas_size=(1800,500),x_lims=nothing, y_lims=nothing,
     sp_markersize=2, vs_title="Visium", sp_title="Cartana", order=true, scale = true)
     vs_count=deepcopy(vs.normCount)
@@ -1353,7 +1369,7 @@ function overlay_visium_cartana_gene(vs::Union{VisiumObject, SlideseqObject}, sp
     MK.current_figure()
 end
 
-function plot_gene_group_spatial(sp_list::Union{Vector{ImagingSpatialObject}, Vector{CartanaObject},Vector{XeniumObject}, Vector{MerfishObject}, Vector{SlideseqObject}, Vector{seqFishObject}, Vector{STARmapObject}}, n_bin, gene_list; group_names::Union{Vector, String, Nothing}=nothing,
+function plot_gene_group_spatial(sp_list::Union{Vector{ImagingSpatialObject}, Vector{CartanaObject},Vector{XeniumObject}, Vector{MerfishObject}, Vector{SlideseqObject}, Vector{seqFishObject}, Vector{STARmapObject}, Vector{StereoSeqObject}}, n_bin, gene_list; group_names::Union{Vector, String, Nothing}=nothing,
     color_range::Vector=["white", "ivory","gold","orange","tomato","red"],legend_min::Union{Float64, Int64}=0, 
     legend_max::Union{Float64, Int64}=1, assay_use = "measured", imp_type = "SpaGE")
     n_obj = length(sp_list)
diff --git a/src/spatial/sp_processing.jl b/src/spatial/sp_processing.jl
index f73cd9f..3181374 100644
--- a/src/spatial/sp_processing.jl
+++ b/src/spatial/sp_processing.jl
@@ -104,12 +104,26 @@ function generate_polygon_counts(sp::AbstractImagingObj)
     join_df = innerjoin(coord_molecules[!,[:gene,:cell]],anno[!,[:number,:cell]], on=:cell)
     join_df.gene = string.(join_df.gene)
     gdf = groupby(join_df, :gene)
-    new_df = DataFrame()
-    for (i, df) in enumerate(gdf)
-        map_dict = countmap(df.number)
-        anno1 = DataFrame(cell=collect(keys(map_dict)),count=collect(values(map_dict)))
-        anno1.gene .= gdf[i].gene[1]
-        new_df = [new_df;anno1]
+    new_df = []
+    p = Progress(length(gdf), dt=0.5, barglyphs=BarGlyphs("[=> ]"), barlen=50, color=:blue)
+    if isa(sp, StereoSeqObject)
+        for df in gdf
+            anno1 = DataFrames.combine(groupby(df, :number), :MIDCount => sum => :count)
+            anno1.gene .= df.gene[1]
+            push!(new_df, anno1)
+            next!(p)
+        end
+        new_df = vcat(new_df...)
+        rename!(new_df, :number => :cell)    
+    else
+        for df in gdf
+            map_dict = countmap(df.number)
+            anno1 = DataFrame(cell=collect(keys(map_dict)), count=collect(values(map_dict)))
+            anno1.gene .= df.gene[1]
+            push!(new_df, anno1)
+            next!(p)
+        end
+        new_df = vcat(new_df...) 
     end
     final_df = unstack(new_df, :cell, :gene, :count)
     final_df .= ifelse.(isequal.(final_df, missing), 0, final_df)
@@ -129,3 +143,63 @@ function generate_polygon_counts(sp::AbstractImagingObj)
     println("Polygon counts was normalized!")
     return sp
 end
+
+function process_scs_file(file_path, add_x, add_y, prefix)
+    df = CSV.File(file_path; delim='\t', header=false) |> DataFrame
+    rename!(df, [:coord, :cell])
+    coords = split.(df[:, :coord], ':')
+    df.x = parse.(Int, getindex.(coords, 1)) .+ add_x
+    df.y = parse.(Int, getindex.(coords, 2)) .+ add_y
+    df.cell = prefix .* "_" .* string.(df.cell)
+    return df
+end
+
+function process_scs_directory(directory)
+    filenames = glob("spot2cell_*.txt", directory)
+    all_dataframes = DataFrame[]
+    p = Progress(length(filenames), desc="Processing Files")
+    for filename in filenames
+        numbers = [parse(Int, m.match) for m in eachmatch(r"\d+", basename(filename))]
+        add_x, add_y = numbers[2], numbers[3]
+        file_path = joinpath(directory, filename)
+        prefix = join(numbers[2:4], ':')
+        df = process_scs_file(file_path, add_x, add_y, prefix)
+        push!(all_dataframes, df)
+        next!(p)
+    end
+    return vcat(all_dataframes...)
+end
+
+function create_stereoseq_scs(scs_results, spot_coord; prefix="sp", min_gene=0, min_cell=0)
+    @info("1. Reading input files...")
+    final_dataframe = process_scs_directory(scs_results)
+    orig_cord = CSV.File(spot_coord; delim='\t') |> DataFrame
+    final_dataframe[!, :spot_loc] = [string(i) * "_" * string(j) for (i, j) in zip(final_dataframe.x, final_dataframe.y)]
+    orig_cord[!, :spot_loc] = [string(i) * "_" * string(j) for (i, j) in zip(orig_cord.x, orig_cord.y)]
+    orig_cord2 = filter(:spot_loc => ∈(Set(final_dataframe.spot_loc)), orig_cord)
+    new_coord = innerjoin(orig_cord2, final_dataframe[:, [:spot_loc, :cell]], on=:spot_loc)
+    gdf = groupby(new_coord, :cell)
+    new_df = []
+    @info("2. Creating count matrix...")
+    p = Progress(length(gdf), dt=0.5, barglyphs=BarGlyphs("[=> ]"), barlen=50, color=:blue)
+    for df in gdf
+        anno1 = DataFrames.combine(groupby(df, :geneID), :MIDCount => sum => :count)
+        anno1.cell .= df.cell[1]
+        push!(new_df, anno1)
+        next!(p)
+    end
+    new_df = vcat(new_df...)
+    @info("3. Constructing StereoSeq object...")
+    final_df = unstack(new_df, :cell, :geneID, :count)
+    final_df .= ifelse.(isequal.(final_df, missing), 0, final_df)
+    cellnames = final_df.cell
+    final_df = mapcols(ByRow(Int64), final_df[!, 2:end])
+    genenames = names(final_df)
+    final_df = convert(SparseMatrixCSC{Int64, Int64},Matrix(final_df)')
+    raw_count = RawCountObject(final_df, cellnames, genenames)
+    cell_coord = DataFrames.combine(groupby(new_coord, :cell), :x => mean => :x, :y => mean => :y)
+    rename!(new_coord,:geneID => :gene)
+    sp = StereoSeqObject(new_coord, cell_coord, raw_count;
+    prefix = prefix, min_gene = min_gene, min_cell = min_cell)
+    return sp
+end
diff --git a/src/spatial/sp_utils.jl b/src/spatial/sp_utils.jl
index 07cb5bd..25cd56e 100644
--- a/src/spatial/sp_utils.jl
+++ b/src/spatial/sp_utils.jl
@@ -41,7 +41,7 @@ function reorder(df::DataFrame,
     return df
 end
 
-function compute_pearson_cor(sp::Union{ImagingSpatialObject, CartanaObject, XeniumObject, MerfishObject, seqFishObject, STARmapObject}, cluster1::Union{Symbol, String}, cluster2::Union{Symbol, String}; color_scheme::String="lightgreyred",reverse_color::Bool=false)
+function compute_pearson_cor(sp::Union{ImagingSpatialObject, CartanaObject, XeniumObject, MerfishObject, seqFishObject, STARmapObject, StereoSeqObject}, cluster1::Union{Symbol, String}, cluster2::Union{Symbol, String}; color_scheme::String="lightgreyred",reverse_color::Bool=false)
     df=sp.spmetaData.cell
     celltypes1=unique(df[!,cluster1])
     celltypes2=unique(df[!,cluster2])
@@ -231,7 +231,7 @@ function slope2deg(slope::Float64)
     return degree
 end
 
-function subset_fov(sp::Union{ImagingSpatialObject, CartanaObject, XeniumObject,MerfishObject, seqFishObject, STARmapObject}, fov::Vector{Int64}, n_fields_x::Int64, n_fields_y::Int64)
+function subset_fov(sp::Union{ImagingSpatialObject, CartanaObject, XeniumObject,MerfishObject, seqFishObject, STARmapObject, StereoSeqObject}, fov::Vector{Int64}, n_fields_x::Int64, n_fields_y::Int64)
     df=sp.spmetaData.cell
     pts, centroids = split_field(df, n_fields_x, n_fields_y)
     pts_sub=pts[fov]
@@ -274,7 +274,7 @@ function compute_new_coord(df, pt, center; span=150)
     return depth, angle
 end
 
-function compute_kidney_coordinates(sp::Union{ImagingSpatialObject, CartanaObject, XeniumObject, MerfishObject, seqFishObject, STARmapObject}, center)
+function compute_kidney_coordinates(sp::Union{ImagingSpatialObject, CartanaObject, XeniumObject, MerfishObject, seqFishObject, STARmapObject, StereoSeqObject}, center)
     df = sp.spmetaData.cell
     kid_depth=[]
     kid_angle=[]
@@ -398,7 +398,7 @@ function visium_unit_radius(spot_num)
     return r_unit
 end
 
-function visium_deconvolution(vs::Union{VisiumObject, SlideseqObject}, sp::Union{ImagingSpatialObject, CartanaObject, XeniumObject, MerfishObject, seqFishObject, STARmapObject}, spot_r::Union{Int64, Float64};
+function visium_deconvolution(vs::Union{VisiumObject, SlideseqObject}, sp::Union{ImagingSpatialObject, CartanaObject, XeniumObject, MerfishObject, seqFishObject, STARmapObject, StereoSeqObject}, spot_r::Union{Int64, Float64};
     vscell_col = "cell" , spcluster_col="celltype", vs_x = "new_x", vs_y = "new_y", 
     sp_x = "new_x", sp_y = "new_y")
     vs_cells = deepcopy(vs.cells)