Mercurial > repos > PredictPDPS / changeset
changeset
Added PET and updated README
Fri, 19 Apr 2024 17:00:37 +0300
- author
- Neil Dizon <neil.dizon@helsinki.fi>
- date
- Fri, 19 Apr 2024 17:00:37 +0300
- changeset 8
- e4ad8f7ce671
- parent 7
- 4bbb246f4cac
- child 9
- 4eaff3502eea
Added PET and updated README
--- a/Manifest.toml Thu Apr 18 11:31:32 2024 +0300 +++ b/Manifest.toml Fri Apr 19 17:00:37 2024 +0300 @@ -13,9 +13,9 @@ [[Adapt]] deps = ["LinearAlgebra", "Requires"] -git-tree-sha1 = "6a55b747d1812e699320963ffde36f1ebdda4099" +git-tree-sha1 = "cde29ddf7e5726c9fb511f340244ea3481267608" uuid = "79e6a3ab-5dfb-504d-930d-738a2a938a0e" -version = "4.0.4" +version = "3.7.2" weakdeps = ["StaticArrays"] [Adapt.extensions] @@ -31,36 +31,20 @@ uuid = "0dad84c5-d112-42e6-8d28-ef12dabb789f" version = "1.1.1" -[[ArrayInterface]] -deps = ["Adapt", "LinearAlgebra", "Requires", "SparseArrays", "SuiteSparse"] -git-tree-sha1 = "c5aeb516a84459e0318a02507d2261edad97eb75" -uuid = "4fba245c-0d91-5ea0-9b3e-6abc04ee57a9" -version = "7.7.1" - - [ArrayInterface.extensions] - ArrayInterfaceBandedMatricesExt = "BandedMatrices" - ArrayInterfaceBlockBandedMatricesExt = "BlockBandedMatrices" - ArrayInterfaceCUDAExt = "CUDA" - ArrayInterfaceGPUArraysCoreExt = "GPUArraysCore" - ArrayInterfaceStaticArraysCoreExt = "StaticArraysCore" - ArrayInterfaceTrackerExt = "Tracker" - - [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" +[[ArnoldiMethod]] +deps = ["LinearAlgebra", "Random", "StaticArrays"] +git-tree-sha1 = "d57bd3762d308bded22c3b82d033bff85f6195c6" +uuid = "ec485272-7323-5ecc-a04f-4719b315124d" +version = "0.4.0" [[Artifacts]] uuid = "56f22d72-fd6d-98f1-02f0-08ddc0907c33" [[AxisAlgorithms]] deps = ["LinearAlgebra", "Random", "SparseArrays", "WoodburyMatrices"] -git-tree-sha1 = "01b8ccb13d68535d73d2b0c23e39bd23155fb712" +git-tree-sha1 = "66771c8d21c8ff5e3a93379480a2307ac36863f7" uuid = "13072b0f-2c55-5437-9ae7-d433b7a33950" -version = "1.1.0" +version = "1.0.1" [[AxisArrays]] deps = ["Dates", "IntervalSets", "IterTools", "RangeArrays"] @@ -76,12 +60,6 @@ uuid = "d1d4a3ce-64b1-5f1a-9ba4-7e7e69966f35" version = "0.1.8" -[[BitTwiddlingConvenienceFunctions]] -deps = ["Static"] -git-tree-sha1 = "0c5f81f47bbbcf4aea7b2959135713459170798b" -uuid = "62783981-4cbd-42fc-bca8-16325de8dc4b" -version = "0.1.5" - [[Bzip2_jll]] deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"] git-tree-sha1 = "9e2a6b69137e6969bab0152632dcb3bc108c8bdd" @@ -93,12 +71,6 @@ uuid = "fa961155-64e5-5f13-b03f-caf6b980ea82" version = "0.5.0" -[[CPUSummary]] -deps = ["CpuId", "IfElse", "PrecompileTools", "Static"] -git-tree-sha1 = "601f7e7b3d36f18790e2caf83a882d88e9b71ff1" -uuid = "2a0fbf3d-bb9c-48f3-b0a9-814d99fd7ab9" -version = "0.2.4" - [[Cairo_jll]] deps = ["Artifacts", "Bzip2_jll", "CompilerSupportLibraries_jll", "Fontconfig_jll", "FreeType2_jll", "Glib_jll", "JLLWrappers", "LZO_jll", "Libdl", "Pixman_jll", "Xorg_libXext_jll", "Xorg_libXrender_jll", "Zlib_jll", "libpng_jll"] git-tree-sha1 = "a4c43f59baa34011e303e76f5c8c91bf58415aaf" @@ -121,11 +93,11 @@ [ChainRulesCore.extensions] ChainRulesCoreSparseArraysExt = "SparseArrays" -[[CloseOpenIntervals]] -deps = ["Static", "StaticArrayInterface"] -git-tree-sha1 = "70232f82ffaab9dc52585e0dd043b5e0c6b714f1" -uuid = "fb6a15b2-703c-40df-9091-08a04967cfa9" -version = "0.1.12" +[[Clustering]] +deps = ["Distances", "LinearAlgebra", "NearestNeighbors", "Printf", "Random", "SparseArrays", "Statistics", "StatsBase"] +git-tree-sha1 = "9ebb045901e9bbf58767a9f34ff89831ed711aae" +uuid = "aaaa29a8-35af-508c-8bc3-b662a17a0fe5" +version = "0.15.7" [[CodecZlib]] deps = ["TranscodingStreams", "Zlib_jll"] @@ -133,6 +105,12 @@ uuid = "944b1d66-785c-5afd-91f1-9de20f533193" version = "0.7.4" +[[ColorSchemes]] +deps = ["ColorTypes", "ColorVectorSpace", "Colors", "FixedPointNumbers", "PrecompileTools", "Random"] +git-tree-sha1 = "67c1f244b991cad9b0aa4b7540fb758c2488b129" +uuid = "35d6a980-a343-548e-a6ea-1d62b119f2f4" +version = "3.24.0" + [[ColorTypes]] deps = ["FixedPointNumbers", "Random"] git-tree-sha1 = "b10d0b65641d57b8b4d5e234446582de5047050d" @@ -177,6 +155,12 @@ uuid = "f0e56b4a-5159-44fe-b623-3e5288b988bb" version = "2.4.1" +[[Conda]] +deps = ["Downloads", "JSON", "VersionParsing"] +git-tree-sha1 = "51cab8e982c5b598eea9c8ceaced4b58d9dd37c9" +uuid = "8f4d0f93-b110-5947-807f-2305c1781a2d" +version = "1.10.0" + [[ConstructionBase]] deps = ["LinearAlgebra"] git-tree-sha1 = "260fd2400ed2dab602a7c15cf10c1933c59930a2" @@ -194,17 +178,16 @@ uuid = "150eb455-5306-5404-9cee-2592286d6298" version = "0.6.3" -[[CpuId]] -deps = ["Markdown"] -git-tree-sha1 = "fcbb72b032692610bfbdb15018ac16a36cf2e406" -uuid = "adafc99b-e345-5852-983c-f28acb93d879" -version = "0.3.1" - [[CustomUnitRanges]] git-tree-sha1 = "1a3f97f907e6dd8983b744d2642651bb162a3f7a" uuid = "dc8bdbbb-1ca9-579f-8c36-e416f6a65cce" version = "1.0.2" +[[DataAPI]] +git-tree-sha1 = "abe83f3a2f1b857aac70ef8b269080af17764bbe" +uuid = "9a962f9c-6df0-11e9-0e5d-c546b8b5ee8a" +version = "1.16.0" + [[DataStructures]] deps = ["Compat", "InteractiveUtils", "OrderedCollections"] git-tree-sha1 = "1d0a14036acb104d9e89698bd408f63ab58cdc82" @@ -374,6 +357,12 @@ uuid = "3b182d85-2403-5c21-9c21-1e1f0cc25472" version = "1.3.14+0" +[[Graphs]] +deps = ["ArnoldiMethod", "Compat", "DataStructures", "Distributed", "Inflate", "LinearAlgebra", "Random", "SharedArrays", "SimpleTraits", "SparseArrays", "Statistics"] +git-tree-sha1 = "3863330da5466410782f2bffc64f3d505a6a8334" +uuid = "86223c79-3864-5bf0-83f7-82e725a168b6" +version = "1.10.0" + [[HTTP]] deps = ["Base64", "CodecZlib", "ConcurrentUtilities", "Dates", "ExceptionUnwrapping", "Logging", "LoggingExtras", "MbedTLS", "NetworkOptions", "OpenSSL", "Random", "SimpleBufferStream", "Sockets", "URIs", "UUIDs"] git-tree-sha1 = "8e59b47b9dc525b70550ca082ce85bcd7f5477cd" @@ -386,17 +375,6 @@ uuid = "2e76f6c2-a576-52d4-95c1-20adfe4de566" version = "2.8.1+1" -[[HostCPUFeatures]] -deps = ["BitTwiddlingConvenienceFunctions", "IfElse", "Libdl", "Static"] -git-tree-sha1 = "eb8fed28f4994600e29beef49744639d985a04b2" -uuid = "3e5b6fbb-0976-4d2c-9146-d79de83f2fb0" -version = "0.1.16" - -[[IfElse]] -git-tree-sha1 = "debdd00ffef04665ccbb3e150747a77560e8fad1" -uuid = "615f187c-cbe4-4ef1-ba3b-2fcf58d6d173" -version = "0.1.1" - [[ImageAxes]] deps = ["AxisArrays", "ImageBase", "ImageCore", "Reexport", "SimpleTraits"] git-tree-sha1 = "2e4520d67b0cef90865b3ef727594d2a58e0e1f8" @@ -458,10 +436,10 @@ version = "0.9.9" [[ImageMorphology]] -deps = ["DataStructures", "ImageCore", "LinearAlgebra", "LoopVectorization", "OffsetArrays", "Requires", "TiledIteration"] -git-tree-sha1 = "6f0a801136cb9c229aebea0df296cdcd471dbcd1" +deps = ["ImageCore", "LinearAlgebra", "Requires", "TiledIteration"] +git-tree-sha1 = "e7c68ab3df4a75511ba33fc5d8d9098007b579a8" uuid = "787d08f9-d448-5407-9aad-5290dd7ab264" -version = "0.4.5" +version = "0.3.2" [[ImageQualityIndexes]] deps = ["ImageContrastAdjustment", "ImageCore", "ImageDistances", "ImageFiltering", "LazyModules", "OffsetArrays", "PrecompileTools", "Statistics"] @@ -469,6 +447,18 @@ uuid = "2996bd0c-7a13-11e9-2da2-2f5ce47296a9" version = "0.3.7" +[[ImageSegmentation]] +deps = ["Clustering", "DataStructures", "Distances", "Graphs", "ImageCore", "ImageFiltering", "ImageMorphology", "LinearAlgebra", "MetaGraphs", "RegionTrees", "SimpleWeightedGraphs", "StaticArrays", "Statistics"] +git-tree-sha1 = "44664eea5408828c03e5addb84fa4f916132fc26" +uuid = "80713f31-8817-5129-9cf8-209ff8fb23e1" +version = "1.8.1" + +[[ImageShow]] +deps = ["Base64", "ColorSchemes", "FileIO", "ImageBase", "ImageCore", "OffsetArrays", "StackViews"] +git-tree-sha1 = "3b5344bcdbdc11ad58f3b1956709b5b9345355de" +uuid = "4e3cecfd-b093-5904-9786-8bbb286a6a31" +version = "0.3.8" + [[ImageTools]] deps = ["AlgTools", "ColorTypes", "FileIO", "GR", "OffsetArrays", "Printf", "Setfield"] path = "../ImageTools" @@ -476,10 +466,16 @@ version = "0.1.0" [[ImageTransformations]] -deps = ["AxisAlgorithms", "CoordinateTransformations", "ImageBase", "ImageCore", "Interpolations", "OffsetArrays", "Rotations", "StaticArrays"] -git-tree-sha1 = "e0884bdf01bbbb111aea77c348368a86fb4b5ab6" +deps = ["AxisAlgorithms", "ColorVectorSpace", "CoordinateTransformations", "ImageBase", "ImageCore", "Interpolations", "OffsetArrays", "Rotations", "StaticArrays"] +git-tree-sha1 = "8717482f4a2108c9358e5c3ca903d3a6113badc9" uuid = "02fcd773-0e25-5acc-982a-7f6622650795" -version = "0.10.1" +version = "0.9.5" + +[[Images]] +deps = ["Base64", "FileIO", "Graphics", "ImageAxes", "ImageBase", "ImageContrastAdjustment", "ImageCore", "ImageDistances", "ImageFiltering", "ImageIO", "ImageMagick", "ImageMetadata", "ImageMorphology", "ImageQualityIndexes", "ImageSegmentation", "ImageShow", "ImageTransformations", "IndirectArrays", "IntegralArrays", "Random", "Reexport", "SparseArrays", "StaticArrays", "Statistics", "StatsBase", "TiledIteration"] +git-tree-sha1 = "5fa9f92e1e2918d9d1243b1131abe623cdf98be7" +uuid = "916415d5-f1e6-5110-898d-aaa5f9f070e0" +version = "0.25.3" [[Imath_jll]] deps = ["Artifacts", "JLLWrappers", "Libdl"] @@ -497,6 +493,12 @@ uuid = "d25df0c9-e2be-5dd7-82c8-3ad0b3e990b9" version = "0.1.4" +[[IntegralArrays]] +deps = ["ColorTypes", "FixedPointNumbers", "IntervalSets"] +git-tree-sha1 = "be8e690c3973443bec584db3346ddc904d4884eb" +uuid = "1d092043-8f09-5a30-832f-7509e371ab51" +version = "0.1.5" + [[IntelOpenMP_jll]] deps = ["Artifacts", "JLLWrappers", "Libdl"] git-tree-sha1 = "5fdf2fe6724d8caabf43b557b84ce53f3b7e2f6b" @@ -509,15 +511,9 @@ [[Interpolations]] deps = ["Adapt", "AxisAlgorithms", "ChainRulesCore", "LinearAlgebra", "OffsetArrays", "Random", "Ratios", "Requires", "SharedArrays", "SparseArrays", "StaticArrays", "WoodburyMatrices"] -git-tree-sha1 = "88a101217d7cb38a7b481ccd50d21876e1d1b0e0" +git-tree-sha1 = "721ec2cf720536ad005cb38f50dbba7b02419a15" uuid = "a98d9a8b-a2ab-59e6-89dd-64a1c18fca59" -version = "0.15.1" - - [Interpolations.extensions] - InterpolationsUnitfulExt = "Unitful" - - [Interpolations.weakdeps] - Unitful = "1986cc42-f94f-5a68-af5c-568840ba703d" +version = "0.14.7" [[IntervalSets]] git-tree-sha1 = "dba9ddf07f77f60450fe5d2e2beb9854d9a49bd0" @@ -544,6 +540,12 @@ uuid = "c8e1da08-722c-5040-9ed9-7db0dc04731e" version = "1.10.0" +[[JLD2]] +deps = ["FileIO", "MacroTools", "Mmap", "OrderedCollections", "Pkg", "PrecompileTools", "Printf", "Reexport", "Requires", "TranscodingStreams", "UUIDs"] +git-tree-sha1 = "5ea6acdd53a51d897672edb694e3cc2912f3f8a7" +uuid = "033835bb-8acc-5ee8-8aae-3f567f8a3819" +version = "0.4.46" + [[JLLWrappers]] deps = ["Artifacts", "Preferences"] git-tree-sha1 = "7e5d6779a1e09a36db2a7b6cff50942a0a7d0fca" @@ -592,11 +594,10 @@ uuid = "dd4b983a-f0e5-5f8d-a1b7-129d4a5fb1ac" version = "2.10.1+0" -[[LayoutPointers]] -deps = ["ArrayInterface", "LinearAlgebra", "ManualMemory", "SIMDTypes", "Static", "StaticArrayInterface"] -git-tree-sha1 = "62edfee3211981241b57ff1cedf4d74d79519277" -uuid = "10f19ff3-798f-405d-979b-55457f8fc047" -version = "0.1.15" +[[LaTeXStrings]] +git-tree-sha1 = "50901ebc375ed41dbf8058da26f9de442febbbec" +uuid = "b964fa9f-0449-5b57-a5c2-d3ea65f4040f" +version = "1.3.1" [[LazyArtifacts]] deps = ["Artifacts", "Pkg"] @@ -712,21 +713,6 @@ uuid = "e6f89c97-d47a-5376-807f-9c37f3926c36" version = "1.0.3" -[[LoopVectorization]] -deps = ["ArrayInterface", "CPUSummary", "CloseOpenIntervals", "DocStringExtensions", "HostCPUFeatures", "IfElse", "LayoutPointers", "LinearAlgebra", "OffsetArrays", "PolyesterWeave", "PrecompileTools", "SIMDTypes", "SLEEFPirates", "Static", "StaticArrayInterface", "ThreadingUtilities", "UnPack", "VectorizationBase"] -git-tree-sha1 = "a13f3be5d84b9c95465d743c82af0b094ef9c2e2" -uuid = "bdcacae8-1622-11e9-2a5c-532679323890" -version = "0.12.169" - - [LoopVectorization.extensions] - ForwardDiffExt = ["ChainRulesCore", "ForwardDiff"] - SpecialFunctionsExt = "SpecialFunctions" - - [LoopVectorization.weakdeps] - ChainRulesCore = "d360d2e6-b24c-11e9-a2a3-2a2ae2dbcce4" - ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210" - SpecialFunctions = "276daf66-3868-5448-9aa4-cd146d93841b" - [[MKL_jll]] deps = ["Artifacts", "IntelOpenMP_jll", "JLLWrappers", "LazyArtifacts", "Libdl"] git-tree-sha1 = "72dc3cf284559eb8f53aa593fe62cb33f83ed0c0" @@ -739,11 +725,6 @@ uuid = "1914dd2f-81c6-5fcd-8719-6d5c9610ff09" version = "0.5.13" -[[ManualMemory]] -git-tree-sha1 = "bcaef4fc7a0cfe2cba636d84cda54b5e4e4ca3cd" -uuid = "d125e4d3-2237-4719-b19c-fa641b8a4667" -version = "0.1.8" - [[MappedArrays]] git-tree-sha1 = "2dab0221fe2b0f2cb6754eaa743cc266339f527e" uuid = "dbb5928d-eab1-5f90-85c2-b9b0edb7c900" @@ -764,6 +745,18 @@ uuid = "c8ffd9c3-330d-5841-b78e-0817d7145fa1" version = "2.28.2+0" +[[MetaGraphs]] +deps = ["Graphs", "JLD2", "Random"] +git-tree-sha1 = "1130dbe1d5276cb656f6e1094ce97466ed700e5a" +uuid = "626554b9-1ddb-594c-aa3c-2596fe9399a5" +version = "0.7.2" + +[[Missings]] +deps = ["DataAPI"] +git-tree-sha1 = "ec4f7fbeab05d7747bdf98eb74d130a2a2ed298d" +uuid = "e1d29d7a-bbdc-5cf2-9ac0-f12de2c33e28" +version = "1.2.0" + [[Mmap]] uuid = "a63ad114-7e13-5084-954f-fe012c677804" @@ -783,6 +776,12 @@ uuid = "77ba4419-2d1f-58cd-9bb1-8ffee604a2e3" version = "1.0.2" +[[NearestNeighbors]] +deps = ["Distances", "StaticArrays"] +git-tree-sha1 = "ded64ff6d4fdd1cb68dfcbb818c69e144a5b2e4c" +uuid = "b8a86587-4115-5ab1-83bc-aa920d37bbce" +version = "0.4.16" + [[Netpbm]] deps = ["FileIO", "ImageCore", "ImageMetadata"] git-tree-sha1 = "d92b107dbb887293622df7697a2223f9f8176fcd" @@ -894,6 +893,12 @@ uuid = "69de0a69-1ddd-5017-9359-2bf0b02dc9f0" version = "2.8.1" +[[PerceptualColourMaps]] +deps = ["Colors", "Images", "Interpolations", "LinearAlgebra", "Printf", "PyPlot", "Statistics", "Test"] +git-tree-sha1 = "cff82e9e59864573aa3f4694cc701e7083becdb8" +uuid = "54e51dfa-9dd7-5231-aa84-a4037b83483a" +version = "0.3.6" + [[Pixman_jll]] deps = ["Artifacts", "CompilerSupportLibraries_jll", "JLLWrappers", "LLVMOpenMP_jll", "Libdl"] git-tree-sha1 = "64779bc4c9784fee475689a1752ef4d5747c5e87" @@ -911,11 +916,11 @@ uuid = "eebad327-c553-4316-9ea0-9fa01ccd7688" version = "0.3.3" -[[PolyesterWeave]] -deps = ["BitTwiddlingConvenienceFunctions", "CPUSummary", "IfElse", "Static", "ThreadingUtilities"] -git-tree-sha1 = "240d7170f5ffdb285f9427b92333c3463bf65bf6" -uuid = "1d0040c9-8b98-4ee7-8388-3f51789ca0ad" -version = "0.2.1" +[[PoissonRandom]] +deps = ["Random"] +git-tree-sha1 = "a0f1159c33f846aa77c3f30ebbc69795e5327152" +uuid = "e409e4f3-bfea-5376-8464-e040bb5c01ab" +version = "0.4.4" [[PrecompileTools]] deps = ["Preferences"] @@ -939,6 +944,18 @@ uuid = "92933f4c-e287-5a05-a399-4b506db050ca" version = "1.10.0" +[[PyCall]] +deps = ["Conda", "Dates", "Libdl", "LinearAlgebra", "MacroTools", "Serialization", "VersionParsing"] +git-tree-sha1 = "9816a3826b0ebf49ab4926e2b18842ad8b5c8f04" +uuid = "438e738f-606a-5dbb-bf0a-cddfbfd45ab0" +version = "1.96.4" + +[[PyPlot]] +deps = ["Colors", "LaTeXStrings", "PyCall", "Sockets", "Test", "VersionParsing"] +git-tree-sha1 = "9220a9dae0369f431168d60adab635f88aca7857" +uuid = "d330b81b-6aea-500a-939a-2ce795aea3ee" +version = "2.11.2" + [[QOI]] deps = ["ColorTypes", "FileIO", "FixedPointNumbers"] git-tree-sha1 = "18e8f4d1426e965c7b532ddd260599e1510d26ce" @@ -997,6 +1014,12 @@ uuid = "189a3867-3050-52da-a836-e630ba90ab69" version = "1.2.2" +[[RegionTrees]] +deps = ["IterTools", "LinearAlgebra", "StaticArrays"] +git-tree-sha1 = "4618ed0da7a251c7f92e869ae1a19c74a7d2a7f9" +uuid = "dee08c22-ab7f-5625-9660-a9af2021b33f" +version = "0.3.2" + [[Requires]] deps = ["UUIDs"] git-tree-sha1 = "838a3a4188e2ded87a4f9f184b4b0d78a1e91cb7" @@ -1019,17 +1042,6 @@ uuid = "ea8e919c-243c-51af-8825-aaa63cd721ce" version = "0.7.0" -[[SIMDTypes]] -git-tree-sha1 = "330289636fb8107c5f32088d2741e9fd7a061a5c" -uuid = "94e857df-77ce-4151-89e5-788b33177be4" -version = "0.1.0" - -[[SLEEFPirates]] -deps = ["IfElse", "Static", "VectorizationBase"] -git-tree-sha1 = "3aac6d68c5e57449f5b9b865c9ba50ac2970c4cf" -uuid = "476501e8-09a2-5ece-8869-fb82de89a1fa" -version = "0.6.42" - [[Serialization]] uuid = "9e88b42a-f829-5b0c-bbe9-9e923198166b" @@ -1054,6 +1066,12 @@ uuid = "699a6c99-e7fa-54fc-8d76-47d257e15c1d" version = "0.9.4" +[[SimpleWeightedGraphs]] +deps = ["Graphs", "LinearAlgebra", "Markdown", "SparseArrays"] +git-tree-sha1 = "4b33e0e081a825dbfaf314decf58fa47e53d6acb" +uuid = "47aef6b3-ad0c-573a-a1e2-d07658019622" +version = "1.4.0" + [[Sixel]] deps = ["Dates", "FileIO", "ImageCore", "IndirectArrays", "OffsetArrays", "REPL", "libsixel_jll"] git-tree-sha1 = "2da10356e31327c7096832eb9cd86307a50b1eb6" @@ -1063,6 +1081,12 @@ [[Sockets]] uuid = "6462fe0b-24de-5631-8697-dd941f90decc" +[[SortingAlgorithms]] +deps = ["DataStructures"] +git-tree-sha1 = "66e0a8e672a0bdfca2c3f5937efb8538b9ddc085" +uuid = "a2af1166-a08f-5f64-846c-94a0d3cef48c" +version = "1.2.1" + [[SparseArrays]] deps = ["Libdl", "LinearAlgebra", "Random", "Serialization", "SuiteSparse_jll"] uuid = "2f01184e-e22b-5df5-ae63-d93ebab69eaf" @@ -1089,23 +1113,6 @@ uuid = "cae243ae-269e-4f55-b966-ac2d0dc13c15" version = "0.1.1" -[[Static]] -deps = ["IfElse"] -git-tree-sha1 = "b366eb1eb68075745777d80861c6706c33f588ae" -uuid = "aedffcd0-7271-4cad-89d0-dc628f76c6d3" -version = "0.8.9" - -[[StaticArrayInterface]] -deps = ["ArrayInterface", "Compat", "IfElse", "LinearAlgebra", "PrecompileTools", "Requires", "SparseArrays", "Static", "SuiteSparse"] -git-tree-sha1 = "5d66818a39bb04bf328e92bc933ec5b4ee88e436" -uuid = "0d7ed370-da01-4f52-bd93-41d350b8b718" -version = "1.5.0" -weakdeps = ["OffsetArrays", "StaticArrays"] - - [StaticArrayInterface.extensions] - StaticArrayInterfaceOffsetArraysExt = "OffsetArrays" - StaticArrayInterfaceStaticArraysExt = "StaticArrays" - [[StaticArrays]] deps = ["LinearAlgebra", "PrecompileTools", "Random", "StaticArraysCore"] git-tree-sha1 = "bf074c045d3d5ffd956fa0a461da38a44685d6b2" @@ -1133,16 +1140,18 @@ uuid = "82ae8749-77ed-4fe6-ae5f-f523153014b0" version = "1.7.0" +[[StatsBase]] +deps = ["DataAPI", "DataStructures", "LinearAlgebra", "LogExpFunctions", "Missings", "Printf", "Random", "SortingAlgorithms", "SparseArrays", "Statistics", "StatsAPI"] +git-tree-sha1 = "5cf7606d6cef84b543b483848d4ae08ad9832b21" +uuid = "2913bbd2-ae8a-5f71-8c99-4fb6c76f3a91" +version = "0.34.3" + [[StringDistances]] deps = ["Distances", "StatsAPI"] git-tree-sha1 = "5b2ca70b099f91e54d98064d5caf5cc9b541ad06" uuid = "88034a9c-02f8-509d-84a9-84ec65e18404" version = "0.11.3" -[[SuiteSparse]] -deps = ["Libdl", "LinearAlgebra", "Serialization", "SparseArrays"] -uuid = "4607b0f0-06f3-5cda-b6b1-a6196a1729e9" - [[SuiteSparse_jll]] deps = ["Artifacts", "Libdl", "Pkg", "libblastrampoline_jll"] uuid = "bea87d4a-7f5b-5778-9afe-8cc45184846c" @@ -1174,12 +1183,6 @@ uuid = "5e47fb64-e119-507b-a336-dd2b206d9990" version = "1.8.0" -[[ThreadingUtilities]] -deps = ["ManualMemory"] -git-tree-sha1 = "eda08f7e9818eb53661b3deb74e3159460dfbc27" -uuid = "8290d209-cae3-49c0-8002-c8c24d57dab5" -version = "0.5.2" - [[TiffImages]] deps = ["ColorTypes", "DataStructures", "DocStringExtensions", "FileIO", "FixedPointNumbers", "IndirectArrays", "Inflate", "Mmap", "OffsetArrays", "PkgVersion", "ProgressMeter", "UUIDs"] git-tree-sha1 = "34cc045dd0aaa59b8bbe86c644679bc57f1d5bd0" @@ -1187,10 +1190,10 @@ version = "0.6.8" [[TiledIteration]] -deps = ["OffsetArrays", "StaticArrayInterface"] -git-tree-sha1 = "1176cc31e867217b06928e2f140c90bd1bc88283" +deps = ["OffsetArrays"] +git-tree-sha1 = "5683455224ba92ef59db72d10690690f4a8dc297" uuid = "06e1c1a7-607b-532d-9fad-de7d9aa2abac" -version = "0.5.0" +version = "0.3.1" [[TranscodingStreams]] git-tree-sha1 = "71509f04d045ec714c4748c785a59045c3736349" @@ -1218,11 +1221,10 @@ [[Unicode]] uuid = "4ec0a83e-493e-50e2-b9ac-8f72acf5a8f5" -[[VectorizationBase]] -deps = ["ArrayInterface", "CPUSummary", "HostCPUFeatures", "IfElse", "LayoutPointers", "Libdl", "LinearAlgebra", "SIMDTypes", "Static", "StaticArrayInterface"] -git-tree-sha1 = "ac377f0a248753a1b1d58bbc92a64f5a726dfb71" -uuid = "3d5dd08c-fd9d-11e8-17fa-ed2836048c2f" -version = "0.21.66" +[[VersionParsing]] +git-tree-sha1 = "58d6e80b4ee071f5efd07fda82cb9fbe17200868" +uuid = "81def892-9a0e-5fdd-b105-ffc91e053289" +version = "1.3.0" [[Vulkan_Loader_jll]] deps = ["Artifacts", "JLLWrappers", "Libdl", "Wayland_jll", "Xorg_libX11_jll", "Xorg_libXrandr_jll", "xkbcommon_jll"] @@ -1244,9 +1246,9 @@ [[WoodburyMatrices]] deps = ["LinearAlgebra", "SparseArrays"] -git-tree-sha1 = "c1a7aa6219628fcd757dede0ca95e245c5cd9511" +git-tree-sha1 = "5f24e158cf4cee437052371455fe361f526da062" uuid = "efce3f68-66dc-5838-9240-27a6d6f5f9b6" -version = "1.0.0" +version = "0.5.6" [[XML2_jll]] deps = ["Artifacts", "JLLWrappers", "Libdl", "Libiconv_jll", "Zlib_jll"]
--- a/Project.toml Thu Apr 18 11:31:32 2024 +0300 +++ b/Project.toml Fri Apr 19 17:00:37 2024 +0300 @@ -4,14 +4,26 @@ [deps] AlgTools = "c46e2e78-5339-41fd-a966-983ff60ab8e7" +ColorSchemes = "35d6a980-a343-548e-a6ea-1d62b119f2f4" ColorTypes = "3da002f7-5984-5a60-b8a6-cbb66c0b333f" +Colors = "5ae59095-9a9b-59fe-a467-6f913c188581" +CoordinateTransformations = "150eb455-5306-5404-9cee-2592286d6298" DelimitedFiles = "8bb1440f-4735-579b-a4ab-409b98df4dab" FileIO = "5789e2e9-d7fb-5bc7-8068-2c6fae9b9549" GR = "28b8d3ca-fb5f-59d9-8090-bfdbd6d07a71" +ImageCore = "a09fc81d-aa75-5fe9-8630-4744c3626534" ImageQualityIndexes = "2996bd0c-7a13-11e9-2da2-2f5ce47296a9" ImageTools = "b548cc0d-4ade-417e-bf62-0e39f9d2eee9" +ImageTransformations = "02fcd773-0e25-5acc-982a-7f6622650795" +Images = "916415d5-f1e6-5110-898d-aaa5f9f070e0" +Interpolations = "a98d9a8b-a2ab-59e6-89dd-64a1c18fca59" +OffsetArrays = "6fe1bfb0-de20-5000-8ca7-80f57d26f881" +PerceptualColourMaps = "54e51dfa-9dd7-5231-aa84-a4037b83483a" +PoissonRandom = "e409e4f3-bfea-5376-8464-e040bb5c01ab" Printf = "de0858da-6303-5e67-8744-51eddeeeb8d7" QuartzImageIO = "dca85d43-d64c-5e67-8c65-017450d5d020" +Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c" +Rotations = "6038ab10-8711-5258-84ad-4b1120ba62dc" Setfield = "efcf1570-3423-57d1-acb7-fd33fddbac46" StableRNGs = "860ef19b-820b-49d6-a774-d7a799459cd3" TestImages = "5e47fb64-e119-507b-a336-dd2b206d9990"
--- a/README.md Thu Apr 18 11:31:32 2024 +0300 +++ b/README.md Fri Apr 19 17:00:37 2024 +0300 @@ -1,13 +1,13 @@ # Julia codes for “Predictive online optimisation with applications to optical flow” -These are the Julia codes for the optical flow experiments of the manuscript _“Predictive online optimisation with applications to optical flow”_ by [Tuomo Valkonen](https://tuomov.iki.fi) ([arXiv:2002.03053](https://arxiv.org/abs/2002.03053)). +This version of Julia codes contains the experiments in the manuscript _"Prediction Techniques for Dynamic Imaging with Online Primal-dual Methods"_ built on top of the experiments for the manuscript _“Predictive online optimisation with applications to optical flow”_ by [Tuomo Valkonen](https://tuomov.iki.fi) ([Journal of Mathematical Imaging and Vision](https://link.springer.com/article/10.1007/s10851-020-01000-4)). ## Prerequisites These codes were written for Julia 1.3. The Julia package prequisites are from November 2019 when our experiments were run, and have not been updated to maintain the same environment we used to do the experiments in the manuscript. You may get Julia from [julialang.org](https://julialang.org/). -## Using +## Usage Navigate your unix shell to the directory containing this `README.md` and then run: @@ -21,16 +21,32 @@ > using PredictPDPS -This may take a while as Julia precompiles the code. Then, to generate all the experiments in the manuscript, run: +This may take a while as Julia precompiles the code. Then, to generate all the experiments in the manuscript _“Predictive online optimisation with applications to optical flow”_, run: > batchrun_article() -This will save the results under `img/`. -To see the experiments running visually, and not save the results, run +To see the experiments running visually, and not save the results, run > demo_known1() -or any of `demo_XY()`, where `X`=1,2,3 and `Y`=`known`,`unknown`. -Further parameters and experiments are available via `run_experiments`. See the source code for details. +or any of `demo_XY()`, where `X`=`known`,`unknown` and `Y`=1,2,3. + +Additionally, to generate all the experiments in the manuscript _"Prediction Techniques for Dynamic Imaging with Online Primal-dual Methods"_, run: + + > batchrun_predictors() + +Both will save the results under `img/`. + +To see the experiments running visually, and not save the results, run + + > demo_denoising1() + +or + + > demo_pet1() + +or any of `demo_denoisingZ()` for image stabilisation experiments and `demo_petZ()` for dynamic PET reconstruction where `Z=1` for Dual Scaling, `Z=2` for Greedy, `Z=3` for No Prediction, `Z=4` for Primal Only, `Z=5` for Proximal, and `Z=6` for Rotation predictors. + +See the source code for more details details. To run the data generation multi-threadeadly parallel to the algorithm, set the `JULIA_NUM_THREADS` environment variable to a number larger than one.
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/src/AlgorithmPET.jl Fri Apr 19 17:00:37 2024 +0300 @@ -0,0 +1,242 @@ +#################################################################### +# Predictive online PDPS for optical flow with known velocity field +#################################################################### + +__precompile__() + +module AlgorithmPET + +identifier = "pet_known_orig" + +using Printf + +using AlgTools.Util +import AlgTools.Iterate +using ImageTools.Gradient +using ImageTools.Translate + +using ..Radon +using ImageTransformations +using Images, CoordinateTransformations, Rotations, OffsetArrays +using ImageCore, Interpolations + +using ..OpticalFlow: ImageSize, + Image, + petpdflow! + +######################### +# Iterate initialisation +######################### + + + +function init_rest(x::Image) + imdim=size(x) + + y = zeros(2, imdim...) + Δx = copy(x) + Δy = copy(y) + x̄ = copy(x) + radonx = copy(x) + + return x, y, Δx, Δy, x̄, radonx +end + +function init_iterates(xinit::Image) + return init_rest(copy(xinit)) +end + +function init_iterates(dim::ImageSize) + return init_rest(zeros(dim...)) +end + +######################### +# PETscan related +######################### +function petvalue(x, b, c) + tmp = similar(b) + radon!(tmp, x) + return sum(@. tmp - b*log(tmp+c)) +end + +function petgrad!(res, x, b, c, S) + tmp = similar(b) + radon!(tmp, x) + @. tmp = S .- b/(tmp+c) + backproject!(res, S.*tmp) +end + +function proj_nonneg!(y) + @inbounds @simd for i=1:length(y) + if y[i] < 0 + y[i] = 0 + end + end + return y +end + + + +############ +# Algorithm +############ + +function step_lengths(params, γ, R_K²) + ρ̃₀, τ₀, σ₀, σ̃₀ = params.ρ̃₀, params.τ₀, params.σ₀, params.σ̃₀ + δ = params.δ + ρ = isdefined(params, :phantom_ρ) ? params.phantom_ρ : params.ρ + Λ = params.Λ + Θ = params.dual_flow ? Λ : 1 + + τ = τ₀/γ + @assert(1+γ*τ ≥ Λ) + σ = σ₀*1/(τ*R_K²) + #σ = σ₀*min(1/(τ*R_K²), 1/max(0, τ*R_K²/((1+γ*τ-Λ)*(1-δ))-ρ)) + q = δ*(1+σ*ρ)/Θ + if 1 ≥ q + σ̃ = σ̃₀*σ/q + #ρ̃ = ρ̃₀*max(0, ((Θ*σ)/(2*δ*σ̃^2*(1+σ*ρ))+1/(2σ)-1/σ̃)) + ρ̃ = max(0, (1-q)/(2*σ)) + else + σ̃ = σ̃₀*σ/(q*(1-√(1-1/q))) + ρ̃ = 0 + end + + #println("Step length parameters: τ=$(τ), σ=$(σ), σ̃=$(σ̃), ρ̃=$(ρ̃)") + + return τ, σ, σ̃, ρ̃ +end + +function solve( :: Type{DisplacementT}; + dim :: ImageSize, + iterate = AlgTools.simple_iterate, + params::NamedTuple) where DisplacementT + + ################################ + # Extract and set up parameters + ################################ + + α, ρ = params.α, params.ρ + R_K² = ∇₂_norm₂₂_est² + γ = 1 + # τ, σ, σ̃, ρ̃ = step_lengths(params, γ, R_K²) + λ = params.λ + ω = 1 + c = params.c*ones(params.radondims...) + + ρ̃₀, τ₀, σ₀, σ̃₀ = params.ρ̃₀, params.τ₀, params.σ₀, params.σ̃₀ + + # Update step length parameters + L = 300.0 + τ = τ₀/L + σ = σ₀*(1-τ₀)/(R_K²*τ) + println("Step length parameters: L=$(round(L, digits=4)), τ=$(round(τ, digits=4)), σ=$(round(σ, digits=4))") + + + + δ = params.δ + ρ = isdefined(params, :phantom_ρ) ? params.phantom_ρ : params.ρ + Λ = params.Λ + Θ = params.dual_flow ? Λ : 1 + + q = δ*(1+σ*ρ)/Θ + if 1 ≥ q + σ̃ = σ̃₀*σ/q + #ρ̃ = ρ̃₀*max(0, ((Θ*σ)/(2*δ*σ̃^2*(1+σ*ρ))+1/(2σ)-1/σ̃)) + ρ̃ = max(0, (1-q)/(2*σ)) + else + σ̃ = σ̃₀*σ/(q*(1-√(1-1/q))) + ρ̃ = 0 + end + + ###################### + # Initialise iterates + ###################### + + x, y, Δx, Δy, x̄, r∇ = init_iterates(dim) + + # L = 1.0 + # oldpetgradx = zeros(size(x)...) + # petgradx = zeros(size(x)) + # oldx = ones(size(x)) + + #################### + # Run the algorithm + #################### + # THIS IS THE step function inside iterate_visualise + v = iterate(params) do verbose :: Function, + b :: Image, # noisy_sinogram + v_known :: DisplacementT, + theta_known :: DisplacementT, + b_true :: Image, + S :: Image + + ################################## + # Update the step length parameter + ################################## + # τ = τ₀/L + # σ = σ₀*(1-τ₀)/(R_K²*τ) + # println("Step length parameters: L=$(round(L, digits=4)), τ=$(round(τ, digits=4)), σ=$(round(σ, digits=4))") + + + ################### + # Prediction steps + ################### + + petpdflow!(x, Δx, y, Δy, v_known, theta_known, params.dual_flow) # Old algorithm + #pdflow!(x, Δx, y, Δy, v_known, theta_known, params.dual_flow, 1e-2,1e-2) # Rotation + #pdflow!(x, Δx, y, Δy, v_known, theta_known, params.dual_flow, 1e-2) # Adhoc + #@. oldx = x + + if params.prox_predict + ∇₂!(Δy, x) + @. y = (y + σ̃*Δy)/(1 + σ̃*(ρ̃+ρ/α)) + #@. cc = y + 1000000*σ̃*Δy + #@. y = (y + σ̃*Δy)/(1 + σ̃*(ρ̃+ρ/α)) + (1 - 1/(1 + ρ̃*σ̃))*cc + proj_norm₂₁ball!(y, α) + end + + + ############ + # PDPS step + ############ + + ∇₂ᵀ!(Δx, y) # primal step: + @. x̄ = x # | save old x for over-relax + petgrad!(r∇, x, b, c, S) # | Calculate gradient of fidelity term + + @. x = x-(τ*λ)*r∇-τ*Δx # | + proj_nonneg!(x) # | non-negativity constaint prox + @. x̄ = (1+ω)*x - ω*x̄ # over-relax: x̄ = 2x-x_old + ∇₂!(Δy, x̄) # dual step: + @. y = y + σ*Δy # | + proj_norm₂₁ball!(y, α) # | prox + + ########################################## + # Compute for the local Lipschitz constant + ########################################## + # petgrad!(petgradx, x, b, c, S) + # petgrad!(oldpetgradx, oldx, b, c, S) + # if norm₂(x-oldx)>1e-12 + # L = max(0.9*norm₂(petgradx - oldpetgradx)/norm₂(x-oldx),L) + # end + + ################################ + # Give function value if needed + ################################ + + v = verbose() do + ∇₂!(Δy, x) + value = λ*petvalue(x, b, c) + params.α*norm₂₁(Δy) + value, x, [NaN, NaN], nothing + end + + v + end + + return x, y, v +end + +end # Module + +
--- a/src/ImGenerate.jl Thu Apr 18 11:31:32 2024 +0300 +++ b/src/ImGenerate.jl Fri Apr 19 17:00:37 2024 +0300 @@ -16,11 +16,21 @@ using ImageTools.Translate using ..OpticalFlow: Image, DisplacementConstant, DisplacementFull +using ..Radon # Added for reproducibility import StableRNGs: StableRNG, Random const rng = StableRNG(9182737465) +# Added for PET +import PoissonRandom: pois_rand +import Random: shuffle +import Images: center, warp +import CoordinateTransformations: recenter +import Rotations: RotMatrix +import Interpolations: Flat + + ############## # Our exports ############## @@ -28,7 +38,9 @@ export ImGen, OnlineData, imgen_square, - imgen_shake + imgen_shake, + PetOnlineData, + imgen_shepplogan_radon ################## # Data structures @@ -50,6 +62,17 @@ v_cumul_true :: DisplacementT end +struct PetOnlineData{DisplacementT} + b_true :: Image + sinogram_true :: Image + sinogram_noisy :: Image + v :: DisplacementT + v_true :: DisplacementT + v_cumul_true :: DisplacementT + theta :: DisplacementT # theta = thetaknown, theta_cumul + S :: Image +end + ################### # Shake generation ################### @@ -83,6 +106,21 @@ pixelwise = (shakefn, sz) -> () -> make_const_u(shakefn(), sz) + +function rotatebytheta(params) + r = params.rotation_factor*(2*rand(rng)-1) + return r +end + +function generate_radonmask(params) + imdim = params.radondims + sino_sparse = params.sino_sparsity + numone = Int64(round(sino_sparse*imdim[1]*imdim[2])) + numzero = imdim[1]*imdim[2]-numone + A = shuffle(rng,reshape([ones(numone); zeros(numzero)],(imdim[1],imdim[2]))) + return A +end + ################ # Moving square ################ @@ -141,9 +179,9 @@ :: Type{DisplacementConstant}, datachannel :: Channel{OnlineData{DisplacementConstant}}, params :: NamedTuple) - + # Restart the seed to enable comparison across predictors - Random.seed!(rng,9182737465) + Random.seed!(rng,9182737465) nextv = shake(params) v_true = nextv() @@ -173,4 +211,72 @@ "$(imname)$(sz[1])x$(sz[2])") end + + +######################################################################## +# PETscan +######################################################################## +function generate_radon(im, sz, + :: Type{DisplacementConstant}, + datachannel :: Channel{PetOnlineData{DisplacementConstant}}, + params :: NamedTuple) + + # Restart the seed to enable comparison across predictors + Random.seed!(rng,9182737465) + + nextv = shake(params) + v_true = nextv() + v_cumul = copy(v_true) + + S_true = generate_radonmask(params) + + theta_true = rotatebytheta(params) + theta_cumul = copy(theta_true) + + while true + # Define the transformation matrix + center_point = (sz[1]/2 + v_true[1], sz[2]/2 + v_true[2]) + tform = recenter(RotMatrix(theta_cumul), center_point) + + # Apply the transformation to the image using warp + b_true = copy(warp(im, tform, axes(im), fillvalue=Flat())) + + v = v_true.*(1.0 .+ params.shake_noise_level.*randn(rng,size(v_true)...)) + theta = theta_true*(1.0 + params.rotation_noise_level.*randn(rng)) + + # Generate the true and noisy sinograms + sinogram_true = zeros(params.radondims...) + sinogram_true .*= params.scale + radon!(sinogram_true, b_true) + sinogram_noisy = copy(sinogram_true) + + for i=1:params.radondims[1], j=1:params.radondims[2] + sinogram_noisy[i, j] += pois_rand(rng,params.noise_level) + end + + # Pass data to iteration routine + data = PetOnlineData{DisplacementConstant}(b_true, sinogram_true, sinogram_noisy, v, v_true, v_cumul, [theta, theta_cumul], S_true) + if !put_unless_closed!(datachannel, data) + return + end + + # Next step shake + v_true = nextv() + v_cumul .+= v_true + # Next theta + theta_true = rotatebytheta(params) + theta_cumul += theta_true + # Next sinogram mask + S_true = generate_radonmask(params) + end +end + + +function imgen_shepplogan_radon(origsize,sz) + im = convert(Array{Float64},TestImages.shepp_logan(origsize, highContrast=true)) + dynrange = maximum(im) + return ImGen(curry(generate_radon, im, sz), sz, 1, dynrange, "shepplogan$(sz[1])x$(sz[2])") +end + + end # Module
--- a/src/OpticalFlow.jl Thu Apr 18 11:31:32 2024 +0300 +++ b/src/OpticalFlow.jl Fri Apr 19 17:00:37 2024 +0300 @@ -11,6 +11,15 @@ import ImageTools.Translate using ImageTools.ImFilter +# using ImageTransformations +# using Images, CoordinateTransformations, Rotations, OffsetArrays +# using Interpolations + +import Images: center, warp +import CoordinateTransformations: recenter +import Rotations: RotMatrix +import Interpolations: Flat + ########## # Exports ########## @@ -31,7 +40,8 @@ Gradient, Displacement, DisplacementFull, DisplacementConstant, HornSchunckData, - filter_hs + filter_hs, + petpdflow! ############################################### # Types (several imported from ImageTools.Translate) @@ -154,6 +164,86 @@ y .= C.*y end + +########################## +# PET +########################## +function petflow_interp!(x::AbstractImage, tmp::AbstractImage, u::DisplacementConstant, theta_known::DisplacementConstant; + threads = false) + tmp .= x + center_point = center(x) .+ u + tform = recenter(RotMatrix(theta_known[1]), center_point) + tmp = warp(x, tform, axes(x), fillvalue=Flat()) + x .= tmp +end + +petflow! = petflow_interp! + +function petpdflow!(x, Δx, y, Δy, u, theta_known, dual_flow; threads=:none) + if dual_flow + @backgroundif (threads==:outer) begin + petflow!(x, Δx, u, theta_known; threads=(threads==:inner)) + end begin + petflow!(@view(y[1, :, :]), @view(Δy[1, :, :]), u, theta_known; threads=(threads==:inner)) + petflow!(@view(y[2, :, :]), @view(Δy[2, :, :]), u, theta_known; threads=(threads==:inner)) + end + else + petflow!(x, Δx, u, theta_known) + end +end + +# Method for greedy predictor +function petpdflow!(x, Δx, y, Δy, u, theta_known, dual_flow, β; threads=:none) + oldx = copy(x) + center_point = center(x) .+ u + tform = recenter(RotMatrix(theta_known[1]), center_point) + Δx = warp(x, tform, axes(x), fillvalue=Flat()) + @. x = Δx + @. Δy = y + if dual_flow + Dxx = copy(Δy) + DΔx = copy(Δy) + ∇₂!(Dxx, x) + ∇₂!(DΔx, oldx) + inds = abs.(Dxx) .≤ β + Dxx[inds] .= 1 + DΔx[inds] .= 1 + y .= y.* DΔx ./ Dxx + end +end + +# Method for affine predictor +function petpdflow!(x, Δx, y, u, theta_known, dual_flow; threads=:none) + oldx = copy(x) + center_point = center(x) .+ u + tform = recenter(RotMatrix(theta_known[1]), center_point) + Δx = warp(x, tform, axes(x), fillvalue=Flat()) + @. x = Δx + C = similar(y) + cc = abs.(x-oldx) + if dual_flow + cm = max(1e-12,maximum(cc)) + c = 1 .* (1 .- cc./ cm) .^(10) + C[1,:,:] .= c + C[2,:,:] .= c + y .= C.*y + end +end + +# Method for rotation prediction (exploiting property of inverse rotation) +function petpdflow!(x, Δx, y, Δy, u, theta_known, dual_flow, β₁, β₂; threads=:none) + if dual_flow + @backgroundif (threads==:outer) begin + petflow!(x, Δx, u, theta_known; threads=(threads==:inner)) + end begin + petflow!(@view(y[1, :, :]), @view(Δy[1, :, :]), u, -theta_known; threads=(threads==:inner)) + petflow!(@view(y[2, :, :]), @view(Δy[2, :, :]), u, -theta_known; threads=(threads==:inner)) + end + else + petflow!(x, Δx, u, theta_known) + end +end + ########################## # Linearised optical flow ##########################
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/src/PET/AlgorithmDualScaling.jl Fri Apr 19 17:00:37 2024 +0300 @@ -0,0 +1,180 @@ +#################################################################### +# Predictive online PDPS for optical flow with known velocity field +#################################################################### + +__precompile__() + +module AlgorithmDualScaling + +identifier = "pdps_known_dualscaling" + +using Printf + +using AlgTools.Util +import AlgTools.Iterate +using ImageTools.Gradient +using ImageTools.Translate + +using ..Radon +using ImageTransformations +using Images, CoordinateTransformations, Rotations, OffsetArrays +using ImageCore, Interpolations + +using ..OpticalFlow: ImageSize, + Image, + petpdflow! + +######################### +# Iterate initialisation +######################### + +function init_rest(x::Image) + imdim=size(x) + y = zeros(2, imdim...) + Δx = copy(x) + Δy = copy(y) + x̄ = copy(x) + radonx = copy(x) + return x, y, Δx, Δy, x̄, radonx +end + +function init_iterates(xinit::Image) + return init_rest(copy(xinit)) +end + +function init_iterates(dim::ImageSize) + return init_rest(zeros(dim...)) +end + +######################### +# PETscan related +######################### +function petvalue(x, b, c) + tmp = similar(b) + radon!(tmp, x) + return sum(@. tmp - b*log(tmp+c)) +end + +function petgrad!(res, x, b, c, S) + tmp = similar(b) + radon!(tmp, x) + @. tmp = S .- b/(tmp+c) + backproject!(res, S.*tmp) +end + +function proj_nonneg!(y) + @inbounds @simd for i=1:length(y) + if y[i] < 0 + y[i] = 0 + end + end + return y +end + +############ +# Algorithm +############ + +function solve( :: Type{DisplacementT}; + dim :: ImageSize, + iterate = AlgTools.simple_iterate, + params::NamedTuple) where DisplacementT + + ################################ + # Extract and set up parameters + ################################ + α, ρ = params.α, params.ρ + R_K² = ∇₂_norm₂₂_est² + γ = 1 + L = params.L + τ₀, σ₀ = params.τ₀, params.σ₀ + τ = τ₀/L + σ = σ₀*(1-τ₀)/(R_K²*τ) + + println("Step length parameters: τ=$(τ), σ=$(σ)") + + λ = params.λ + c = params.c*ones(params.radondims...) + + + ###################### + # Initialise iterates + ###################### + + x, y, Δx, Δy, x̄, r∇ = init_iterates(dim) + + if params.L_experiment + oldpetgradx = zeros(size(x)...) + petgradx = zeros(size(x)) + oldx = ones(size(x)) + end + + #################### + # Run the algorithm + #################### + + v = iterate(params) do verbose :: Function, + b :: Image, # noisy_sinogram + v_known :: DisplacementT, + theta_known :: DisplacementT, + b_true :: Image, + S :: Image + + ################### + # Prediction steps + ################### + + petpdflow!(x, Δx, y, v_known, theta_known, true) + + if params.L_experiment + @. oldx = x + end + + ############ + # PDPS step + ############ + + ∇₂ᵀ!(Δx, y) # primal step: + @. x̄ = x # | save old x for over-relax + petgrad!(r∇, x, b, c, S) # | Calculate gradient of fidelity term + + @. x = x-(τ*λ)*r∇-τ*Δx # | + proj_nonneg!(x) # | non-negativity constaint prox + @. x̄ = 2x - x̄ # over-relax: x̄ = 2x-x_old + ∇₂!(Δy, x̄) # dual step: + @. y = y + σ*Δy # | + proj_norm₂₁ball!(y, α) # | prox + + ##################### + # L update if needed + ##################### + if params.L_experiment + petgrad!(petgradx, x, b, c, S) + petgrad!(oldpetgradx, oldx, b, c, S) + if norm₂(x-oldx)>1e-12 + L = max(0.9*norm₂(petgradx - oldpetgradx)/norm₂(x-oldx),L) + end + println("Step length parameters: L=$(L)") + τ = τ₀/L + σ = σ₀*(1-τ₀)/(R_K²*τ) + end + + ################################ + # Give function value if needed + ################################ + + v = verbose() do + ∇₂!(Δy, x) + value = λ*petvalue(x, b, c) + params.α*norm₂₁(Δy) + value, x, [NaN, NaN], nothing + end + + v + end + + return x, y, v +end + +end # Module + +
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/src/PET/AlgorithmGreedy.jl Fri Apr 19 17:00:37 2024 +0300 @@ -0,0 +1,180 @@ +#################################################################### +# Predictive online PDPS for optical flow with known velocity field +#################################################################### + +__precompile__() + +module AlgorithmGreedy + +identifier = "pdps_known_greedy" + +using Printf + +using AlgTools.Util +import AlgTools.Iterate +using ImageTools.Gradient +using ImageTools.Translate + +using ..Radon +using ImageTransformations +using Images, CoordinateTransformations, Rotations, OffsetArrays +using ImageCore, Interpolations + +using ..OpticalFlow: ImageSize, + Image, + petpdflow! + +######################### +# Iterate initialisation +######################### + +function init_rest(x::Image) + imdim=size(x) + y = zeros(2, imdim...) + Δx = copy(x) + Δy = copy(y) + x̄ = copy(x) + radonx = copy(x) + return x, y, Δx, Δy, x̄, radonx +end + +function init_iterates(xinit::Image) + return init_rest(copy(xinit)) +end + +function init_iterates(dim::ImageSize) + return init_rest(zeros(dim...)) +end + +######################### +# PETscan related +######################### +function petvalue(x, b, c) + tmp = similar(b) + radon!(tmp, x) + return sum(@. tmp - b*log(tmp+c)) +end + +function petgrad!(res, x, b, c, S) + tmp = similar(b) + radon!(tmp, x) + @. tmp = S .- b/(tmp+c) + backproject!(res, S.*tmp) +end + +function proj_nonneg!(y) + @inbounds @simd for i=1:length(y) + if y[i] < 0 + y[i] = 0 + end + end + return y +end + +############ +# Algorithm +############ + +function solve( :: Type{DisplacementT}; + dim :: ImageSize, + iterate = AlgTools.simple_iterate, + params::NamedTuple) where DisplacementT + + ################################ + # Extract and set up parameters + ################################ + α, ρ = params.α, params.ρ + R_K² = ∇₂_norm₂₂_est² + γ = 1 + L = params.L + τ₀, σ₀ = params.τ₀, params.σ₀ + τ = τ₀/L + σ = σ₀*(1-τ₀)/(R_K²*τ) + + println("Step length parameters: τ=$(τ), σ=$(σ)") + + λ = params.λ + c = params.c*ones(params.radondims...) + + + ###################### + # Initialise iterates + ###################### + + x, y, Δx, Δy, x̄, r∇ = init_iterates(dim) + + if params.L_experiment + oldpetgradx = zeros(size(x)...) + petgradx = zeros(size(x)) + oldx = ones(size(x)) + end + + #################### + # Run the algorithm + #################### + + v = iterate(params) do verbose :: Function, + b :: Image, # noisy_sinogram + v_known :: DisplacementT, + theta_known :: DisplacementT, + b_true :: Image, + S :: Image + + ################### + # Prediction steps + ################### + + petpdflow!(x, Δx, y, Δy, v_known, theta_known, true, 1e-2) + + if params.L_experiment + @. oldx = x + end + + ############ + # PDPS step + ############ + + ∇₂ᵀ!(Δx, y) # primal step: + @. x̄ = x # | save old x for over-relax + petgrad!(r∇, x, b, c, S) # | Calculate gradient of fidelity term + + @. x = x-(τ*λ)*r∇-τ*Δx # | + proj_nonneg!(x) # | non-negativity constaint prox + @. x̄ = 2x - x̄ # over-relax: x̄ = 2x-x_old + ∇₂!(Δy, x̄) # dual step: + @. y = y + σ*Δy # | + proj_norm₂₁ball!(y, α) # | prox + + ##################### + # L update if needed + ##################### + if params.L_experiment + petgrad!(petgradx, x, b, c, S) + petgrad!(oldpetgradx, oldx, b, c, S) + if norm₂(x-oldx)>1e-12 + L = max(0.9*norm₂(petgradx - oldpetgradx)/norm₂(x-oldx),L) + end + println("Step length parameters: L=$(L)") + τ = τ₀/L + σ = σ₀*(1-τ₀)/(R_K²*τ) + end + + ################################ + # Give function value if needed + ################################ + + v = verbose() do + ∇₂!(Δy, x) + value = λ*petvalue(x, b, c) + params.α*norm₂₁(Δy) + value, x, [NaN, NaN], nothing + end + + v + end + + return x, y, v +end + +end # Module + +
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/src/PET/AlgorithmNoPrediction.jl Fri Apr 19 17:00:37 2024 +0300 @@ -0,0 +1,180 @@ +#################################################################### +# Predictive online PDPS for optical flow with known velocity field +#################################################################### + +__precompile__() + +module AlgorithmNoPrediction + +identifier = "pdps_known_noprediction" + +using Printf + +using AlgTools.Util +import AlgTools.Iterate +using ImageTools.Gradient +using ImageTools.Translate + +using ..Radon +using ImageTransformations +using Images, CoordinateTransformations, Rotations, OffsetArrays +using ImageCore, Interpolations + +using ..OpticalFlow: ImageSize, + Image, + petpdflow! + +######################### +# Iterate initialisation +######################### + +function init_rest(x::Image) + imdim=size(x) + y = zeros(2, imdim...) + Δx = copy(x) + Δy = copy(y) + x̄ = copy(x) + radonx = copy(x) + return x, y, Δx, Δy, x̄, radonx +end + +function init_iterates(xinit::Image) + return init_rest(copy(xinit)) +end + +function init_iterates(dim::ImageSize) + return init_rest(zeros(dim...)) +end + +######################### +# PETscan related +######################### +function petvalue(x, b, c) + tmp = similar(b) + radon!(tmp, x) + return sum(@. tmp - b*log(tmp+c)) +end + +function petgrad!(res, x, b, c, S) + tmp = similar(b) + radon!(tmp, x) + @. tmp = S .- b/(tmp+c) + backproject!(res, S.*tmp) +end + +function proj_nonneg!(y) + @inbounds @simd for i=1:length(y) + if y[i] < 0 + y[i] = 0 + end + end + return y +end + +############ +# Algorithm +############ + +function solve( :: Type{DisplacementT}; + dim :: ImageSize, + iterate = AlgTools.simple_iterate, + params::NamedTuple) where DisplacementT + + ################################ + # Extract and set up parameters + ################################ + α, ρ = params.α, params.ρ + R_K² = ∇₂_norm₂₂_est² + γ = 1 + L = params.L + τ₀, σ₀ = params.τ₀, params.σ₀ + τ = τ₀/L + σ = σ₀*(1-τ₀)/(R_K²*τ) + + println("Step length parameters: τ=$(τ), σ=$(σ)") + + λ = params.λ + c = params.c*ones(params.radondims...) + + + ###################### + # Initialise iterates + ###################### + + x, y, Δx, Δy, x̄, r∇ = init_iterates(dim) + + if params.L_experiment + oldpetgradx = zeros(size(x)...) + petgradx = zeros(size(x)) + oldx = ones(size(x)) + end + + #################### + # Run the algorithm + #################### + + v = iterate(params) do verbose :: Function, + b :: Image, # noisy_sinogram + v_known :: DisplacementT, + theta_known :: DisplacementT, + b_true :: Image, + S :: Image + + ################### + # Prediction steps + ################### + + # No prediction + + if params.L_experiment + @. oldx = x + end + + ############ + # PDPS step + ############ + + ∇₂ᵀ!(Δx, y) # primal step: + @. x̄ = x # | save old x for over-relax + petgrad!(r∇, x, b, c, S) # | Calculate gradient of fidelity term + + @. x = x-(τ*λ)*r∇-τ*Δx # | + proj_nonneg!(x) # | non-negativity constaint prox + @. x̄ = 2x - x̄ # over-relax: x̄ = 2x-x_old + ∇₂!(Δy, x̄) # dual step: + @. y = y + σ*Δy # | + proj_norm₂₁ball!(y, α) # | prox + + ##################### + # L update if needed + ##################### + if params.L_experiment + petgrad!(petgradx, x, b, c, S) + petgrad!(oldpetgradx, oldx, b, c, S) + if norm₂(x-oldx)>1e-12 + L = max(0.9*norm₂(petgradx - oldpetgradx)/norm₂(x-oldx),L) + end + println("Step length parameters: L=$(L)") + τ = τ₀/L + σ = σ₀*(1-τ₀)/(R_K²*τ) + end + + ################################ + # Give function value if needed + ################################ + + v = verbose() do + ∇₂!(Δy, x) + value = λ*petvalue(x, b, c) + params.α*norm₂₁(Δy) + value, x, [NaN, NaN], nothing + end + + v + end + + return x, y, v +end + +end # Module + +
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/src/PET/AlgorithmPrimalOnly.jl Fri Apr 19 17:00:37 2024 +0300 @@ -0,0 +1,180 @@ +#################################################################### +# Predictive online PDPS for optical flow with known velocity field +#################################################################### + +__precompile__() + +module AlgorithmPrimalOnly + +identifier = "pdps_known_primalonly" + +using Printf + +using AlgTools.Util +import AlgTools.Iterate +using ImageTools.Gradient +using ImageTools.Translate + +using ..Radon +using ImageTransformations +using Images, CoordinateTransformations, Rotations, OffsetArrays +using ImageCore, Interpolations + +using ..OpticalFlow: ImageSize, + Image, + petpdflow! + +######################### +# Iterate initialisation +######################### + +function init_rest(x::Image) + imdim=size(x) + y = zeros(2, imdim...) + Δx = copy(x) + Δy = copy(y) + x̄ = copy(x) + radonx = copy(x) + return x, y, Δx, Δy, x̄, radonx +end + +function init_iterates(xinit::Image) + return init_rest(copy(xinit)) +end + +function init_iterates(dim::ImageSize) + return init_rest(zeros(dim...)) +end + +######################### +# PETscan related +######################### +function petvalue(x, b, c) + tmp = similar(b) + radon!(tmp, x) + return sum(@. tmp - b*log(tmp+c)) +end + +function petgrad!(res, x, b, c, S) + tmp = similar(b) + radon!(tmp, x) + @. tmp = S .- b/(tmp+c) + backproject!(res, S.*tmp) +end + +function proj_nonneg!(y) + @inbounds @simd for i=1:length(y) + if y[i] < 0 + y[i] = 0 + end + end + return y +end + +############ +# Algorithm +############ + +function solve( :: Type{DisplacementT}; + dim :: ImageSize, + iterate = AlgTools.simple_iterate, + params::NamedTuple) where DisplacementT + + ################################ + # Extract and set up parameters + ################################ + α, ρ = params.α, params.ρ + R_K² = ∇₂_norm₂₂_est² + γ = 1 + L = params.L + τ₀, σ₀ = params.τ₀, params.σ₀ + τ = τ₀/L + σ = σ₀*(1-τ₀)/(R_K²*τ) + + println("Step length parameters: τ=$(τ), σ=$(σ)") + + λ = params.λ + c = params.c*ones(params.radondims...) + + + ###################### + # Initialise iterates + ###################### + + x, y, Δx, Δy, x̄, r∇ = init_iterates(dim) + + if params.L_experiment + oldpetgradx = zeros(size(x)...) + petgradx = zeros(size(x)) + oldx = ones(size(x)) + end + + #################### + # Run the algorithm + #################### + + v = iterate(params) do verbose :: Function, + b :: Image, # noisy_sinogram + v_known :: DisplacementT, + theta_known :: DisplacementT, + b_true :: Image, + S :: Image + + ################### + # Prediction steps + ################### + + petpdflow!(x, Δx, y, Δy, v_known, theta_known, false) + + if params.L_experiment + @. oldx = x + end + + ############ + # PDPS step + ############ + + ∇₂ᵀ!(Δx, y) # primal step: + @. x̄ = x # | save old x for over-relax + petgrad!(r∇, x, b, c, S) # | Calculate gradient of fidelity term + + @. x = x-(τ*λ)*r∇-τ*Δx # | + proj_nonneg!(x) # | non-negativity constaint prox + @. x̄ = 2x - x̄ # over-relax: x̄ = 2x-x_old + ∇₂!(Δy, x̄) # dual step: + @. y = y + σ*Δy # | + proj_norm₂₁ball!(y, α) # | prox + + ##################### + # L update if needed + ##################### + if params.L_experiment + petgrad!(petgradx, x, b, c, S) + petgrad!(oldpetgradx, oldx, b, c, S) + if norm₂(x-oldx)>1e-12 + L = max(0.9*norm₂(petgradx - oldpetgradx)/norm₂(x-oldx),L) + end + println("Step length parameters: L=$(L)") + τ = τ₀/L + σ = σ₀*(1-τ₀)/(R_K²*τ) + end + + ################################ + # Give function value if needed + ################################ + + v = verbose() do + ∇₂!(Δy, x) + value = λ*petvalue(x, b, c) + params.α*norm₂₁(Δy) + value, x, [NaN, NaN], nothing + end + + v + end + + return x, y, v +end + +end # Module + +
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/src/PET/AlgorithmProximal.jl Fri Apr 19 17:00:37 2024 +0300 @@ -0,0 +1,208 @@ +#################################################################### +# Predictive online PDPS for optical flow with known velocity field +#################################################################### + +__precompile__() + +module AlgorithmProximal + +identifier = "pdps_known_proximal" + +using Printf + +using AlgTools.Util +import AlgTools.Iterate +using ImageTools.Gradient +using ImageTools.Translate + +using ..Radon +using ImageTransformations +using Images, CoordinateTransformations, Rotations, OffsetArrays +using ImageCore, Interpolations + +using ..OpticalFlow: ImageSize, + Image, + petpdflow! + +######################### +# Iterate initialisation +######################### + +function init_rest(x::Image) + imdim=size(x) + y = zeros(2, imdim...) + Δx = copy(x) + Δy = copy(y) + x̄ = copy(x) + radonx = copy(x) + return x, y, Δx, Δy, x̄, radonx +end + +function init_iterates(xinit::Image) + return init_rest(copy(xinit)) +end + +function init_iterates(dim::ImageSize) + return init_rest(zeros(dim...)) +end + +######################### +# PETscan related +######################### +function petvalue(x, b, c) + tmp = similar(b) + radon!(tmp, x) + return sum(@. tmp - b*log(tmp+c)) +end + +function petgrad!(res, x, b, c, S) + tmp = similar(b) + radon!(tmp, x) + @. tmp = S .- b/(tmp+c) + backproject!(res, S.*tmp) +end + +function proj_nonneg!(y) + @inbounds @simd for i=1:length(y) + if y[i] < 0 + y[i] = 0 + end + end + return y +end + +############ +# Algorithm +############ + +function step_lengths(params, γ, R_K², L) + ρ̃₀, τ₀, σ₀, σ̃₀ = params.ρ̃₀, params.τ₀, params.σ₀, params.σ̃₀ + δ = params.δ + ρ = isdefined(params, :phantom_ρ) ? params.phantom_ρ : params.ρ + Λ = params.Λ + Θ = params.dual_flow ? Λ : 1 + + τ = τ₀/L + @assert(1+γ*τ ≥ Λ) + σ = σ₀*min((1-τ₀)/(τ*R_K²), 1/max(0, τ*R_K²/((1+γ*τ-Λ)*(1-δ))-ρ)) + q = δ*(1+σ*ρ)/Θ + if 1 ≥ q + σ̃ = σ̃₀*σ/q + #ρ̃ = ρ̃₀*max(0, ((Θ*σ)/(2*δ*σ̃^2*(1+σ*ρ))+1/(2σ)-1/σ̃)) + ρ̃ = max(0, (1-q)/(2*σ)) + else + σ̃ = σ̃₀*σ/(q*(1-√(1-1/q))) + ρ̃ = 0 + end + + return τ, σ, σ̃, ρ̃ +end + +function solve( :: Type{DisplacementT}; + dim :: ImageSize, + iterate = AlgTools.simple_iterate, + params::NamedTuple) where DisplacementT + + ################################ + # Extract and set up parameters + ################################ + α, ρ = params.α, params.ρ + R_K² = ∇₂_norm₂₂_est² + γ = 1 + L = params.L + τ, σ, σ̃, ρ̃ = step_lengths(params, γ, R_K², L) + println("Step length parameters: τ=$(τ), σ=$(σ), σ̃=$(σ̃), ρ̃=$(ρ̃)") + + λ = params.λ + c = params.c*ones(params.radondims...) + + + + ###################### + # Initialise iterates + ###################### + + x, y, Δx, Δy, x̄, r∇ = init_iterates(dim) + + if params.L_experiment + oldpetgradx = zeros(size(x)...) + petgradx = zeros(size(x)) + oldx = ones(size(x)) + end + + #################### + # Run the algorithm + #################### + + v = iterate(params) do verbose :: Function, + b :: Image, # noisy_sinogram + v_known :: DisplacementT, + theta_known :: DisplacementT, + b_true :: Image, + S :: Image + + ################### + # Prediction steps + ################### + + petpdflow!(x, Δx, y, Δy, v_known, theta_known, params.dual_flow) # Usual flow + + if params.L_experiment + @. oldx = x + end + + if params.prox_predict + ∇₂!(Δy, x) + @. y = (y + σ̃*Δy)/(1 + σ̃*(ρ̃+ρ/α)) + #@. cc = y + 1000000*σ̃*Δy + #@. y = (y + σ̃*Δy)/(1 + σ̃*(ρ̃+ρ/α)) + (1 - 1/(1 + ρ̃*σ̃))*cc + proj_norm₂₁ball!(y, α) + end + + ############ + # PDPS step + ############ + + ∇₂ᵀ!(Δx, y) # primal step: + @. x̄ = x # | save old x for over-relax + petgrad!(r∇, x, b, c, S) # | Calculate gradient of fidelity term + + @. x = x-(τ*λ)*r∇-τ*Δx # | + proj_nonneg!(x) # | non-negativity constaint prox + @. x̄ = 2*x - x̄ # over-relax: x̄ = 2x-x_old + ∇₂!(Δy, x̄) # dual step: + @. y = y + σ*Δy # | + proj_norm₂₁ball!(y, α) # | prox + + ##################### + # L update if needed + ##################### + if params.L_experiment + petgrad!(petgradx, x, b, c, S) + petgrad!(oldpetgradx, oldx, b, c, S) + if norm₂(x-oldx)>1e-12 + L = max(0.9*norm₂(petgradx - oldpetgradx)/norm₂(x-oldx),L) + end + #println("Step length parameters: L=$(L)") + τ, σ, σ̃, ρ̃ = step_lengths(params, γ, R_K², L) + end + + ################################ + # Give function value if needed + ################################ + + v = verbose() do + ∇₂!(Δy, x) + value = λ*petvalue(x, b, c) + params.α*norm₂₁(Δy) + value, x, [NaN, NaN], nothing + end + + v + end + + return x, y, v +end + +end # Module + +
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/src/PET/AlgorithmRotation.jl Fri Apr 19 17:00:37 2024 +0300 @@ -0,0 +1,180 @@ +#################################################################### +# Predictive online PDPS for optical flow with known velocity field +#################################################################### + +__precompile__() + +module AlgorithmRotation + +identifier = "pdps_known_rotation" + +using Printf + +using AlgTools.Util +import AlgTools.Iterate +using ImageTools.Gradient +using ImageTools.Translate + +using ..Radon +using ImageTransformations +using Images, CoordinateTransformations, Rotations, OffsetArrays +using ImageCore, Interpolations + +using ..OpticalFlow: ImageSize, + Image, + petpdflow! + +######################### +# Iterate initialisation +######################### + +function init_rest(x::Image) + imdim=size(x) + y = zeros(2, imdim...) + Δx = copy(x) + Δy = copy(y) + x̄ = copy(x) + radonx = copy(x) + return x, y, Δx, Δy, x̄, radonx +end + +function init_iterates(xinit::Image) + return init_rest(copy(xinit)) +end + +function init_iterates(dim::ImageSize) + return init_rest(zeros(dim...)) +end + +######################### +# PETscan related +######################### +function petvalue(x, b, c) + tmp = similar(b) + radon!(tmp, x) + return sum(@. tmp - b*log(tmp+c)) +end + +function petgrad!(res, x, b, c, S) + tmp = similar(b) + radon!(tmp, x) + @. tmp = S .- b/(tmp+c) + backproject!(res, S.*tmp) +end + +function proj_nonneg!(y) + @inbounds @simd for i=1:length(y) + if y[i] < 0 + y[i] = 0 + end + end + return y +end + +############ +# Algorithm +############ + +function solve( :: Type{DisplacementT}; + dim :: ImageSize, + iterate = AlgTools.simple_iterate, + params::NamedTuple) where DisplacementT + + ################################ + # Extract and set up parameters + ################################ + α, ρ = params.α, params.ρ + R_K² = ∇₂_norm₂₂_est² + γ = 1 + L = params.L + τ₀, σ₀ = params.τ₀, params.σ₀ + τ = τ₀/L + σ = σ₀*(1-τ₀)/(R_K²*τ) + + println("Step length parameters: τ=$(τ), σ=$(σ)") + + λ = params.λ + c = params.c*ones(params.radondims...) + + + ###################### + # Initialise iterates + ###################### + + x, y, Δx, Δy, x̄, r∇ = init_iterates(dim) + + if params.L_experiment + oldpetgradx = zeros(size(x)...) + petgradx = zeros(size(x)) + oldx = ones(size(x)) + end + + #################### + # Run the algorithm + #################### + + v = iterate(params) do verbose :: Function, + b :: Image, # noisy_sinogram + v_known :: DisplacementT, + theta_known :: DisplacementT, + b_true :: Image, + S :: Image + + ################### + # Prediction steps + ################### + + petpdflow!(x, Δx, y, Δy, v_known, theta_known, params.dual_flow, 0,0) + + if params.L_experiment + @. oldx = x + end + + ############ + # PDPS step + ############ + + ∇₂ᵀ!(Δx, y) # primal step: + @. x̄ = x # | save old x for over-relax + petgrad!(r∇, x, b, c, S) # | Calculate gradient of fidelity term + + @. x = x-(τ*λ)*r∇-τ*Δx # | + proj_nonneg!(x) # | non-negativity constaint prox + @. x̄ = 2x - x̄ # over-relax: x̄ = 2x-x_old + ∇₂!(Δy, x̄) # dual step: + @. y = y + σ*Δy # | + proj_norm₂₁ball!(y, α) # | prox + + ##################### + # L update if needed + ##################### + if params.L_experiment + petgrad!(petgradx, x, b, c, S) + petgrad!(oldpetgradx, oldx, b, c, S) + if norm₂(x-oldx)>1e-12 + L = max(0.9*norm₂(petgradx - oldpetgradx)/norm₂(x-oldx),L) + end + println("Step length parameters: L=$(L)") + τ = τ₀/L + σ = σ₀*(1-τ₀)/(R_K²*τ) + end + + ################################ + # Give function value if needed + ################################ + + v = verbose() do + ∇₂!(Δy, x) + value = λ*petvalue(x, b, c) + params.α*norm₂₁(Δy) + value, x, [NaN, NaN], nothing + end + + v + end + + return x, y, v +end + +end # Module + +
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/src/PET/ImGenerate.jl Fri Apr 19 17:00:37 2024 +0300 @@ -0,0 +1,282 @@ +################### +# Image generation +################### + +module ImGenerate + +using ColorTypes: Gray +import TestImages +# We don't really *directly* depend on QuartzImageIO. The import here is +# merely a workaround to suppress warnings when loading TestImages. +# Something is broken in those packages. +import QuartzImageIO + +using AlgTools.Util +using AlgTools.Comms +using ImageTools.Translate + +using ..OpticalFlow: Image, DisplacementConstant, DisplacementFull +using ..Radon + +# Added for reproducibility +import StableRNGs: StableRNG, Random +const rng = StableRNG(9182737465) + +# Added for PET +import PoissonRandom: pois_rand +import Random: shuffle +import Images: center, warp +import CoordinateTransformations: recenter +import Rotations: RotMatrix +import Interpolations: Flat + + +############## +# Our exports +############## + +export ImGen, + OnlineData, + imgen_square, + imgen_shake, + PetOnlineData, + imgen_shepplogan_radon + +################## +# Data structures +################## + +struct ImGen + f :: Function + dim :: Tuple{Int64,Int64} + Λ :: Float64 + dynrange :: Float64 + name :: String +end + +struct OnlineData{DisplacementT} + b_true :: Image + b_noisy :: Image + v :: DisplacementT + v_true :: DisplacementT + v_cumul_true :: DisplacementT +end + +struct PetOnlineData{DisplacementT} + b_true :: Image + sinogram_true :: Image + sinogram_noisy :: Image + v :: DisplacementT + v_true :: DisplacementT + v_cumul_true :: DisplacementT + theta :: DisplacementT # theta = thetaknown, theta_cumul + S :: Image +end + +################### +# Shake generation +################### + +function make_const_v(displ, sz) + v = zeros(2, sz...) + v[1, :, :] .= displ[1] + v[2, :, :] .= displ[2] + return v +end + +function shake(params) + if !haskey(params, :shaketype) || params.shaketype == :gaussian + return () -> params.shake.*randn(rng,2) + elseif params.shaketype == :disk + return () -> begin + θ = 2π*rand(rng,Float64) + r = params.shake*√(rand(rng,Float64)) + return [r*cos(θ), r*sin(θ)] + end + elseif params.shaketype == :circle + return () -> begin + θ = 2π*rand(rng,Float64) + r = params.shake + return [r*cos(θ), r*sin(θ)] + end + else + error("Unknown shaketype $(params.shaketype)") + end +end + +pixelwise = (shakefn, sz) -> () -> make_const_u(shakefn(), sz) + + +function rotatebytheta(params) + r = params.rotation_factor*(2*rand(rng)-1) + return r +end + +function generate_radonmask(params) + imdim = params.radondims + sino_sparse = params.sino_sparsity + numone = Int64(round(sino_sparse*imdim[1]*imdim[2])) + numzero = imdim[1]*imdim[2]-numone + A = shuffle(rng,reshape([ones(numone); zeros(numzero)],(imdim[1],imdim[2]))) + return A +end + +################ +# Moving square +################ + +function generate_square(sz, + :: Type{DisplacementT}, + datachannel :: Channel{OnlineData{DisplacementT}}, + params) where DisplacementT + + if false + v₀ = make_const_v(0.1.*(-1, 1), sz) + nextv = () -> v₀ + elseif DisplacementT == DisplacementFull + nextv = pixelwise(shake(params), sz) + elseif DisplacementT == DisplacementConstant + nextv = shake(params) + else + @error "Invalid DisplacementT" + end + + # Constant linear displacement everywhere has Jacobian determinant one + # (modulo the boundaries which we ignore here) + m = round(Int, sz[1]/5) + b_orig = zeros(sz...) + b_orig[sz[1].-(2*m:3*m), 2*m:3*m] .= 1 + + v_true = nextv() + v_cumul = copy(v_true) + + while true + # Flow original data and add noise + b_true = zeros(sz...) + translate_image!(b_true, b_orig, v_cumul; threads=true) + b = b_true .+ params.noise_level.*randn(rng,sz...) + v = v_true.*(1.0 .+ params.shake_noise_level.*randn(rng,size(v_true)...)) + # Pass true data to iteration routine + data = OnlineData{DisplacementT}(b_true, b, v, v_true, v_cumul) + if !put_unless_closed!(datachannel, data) + return + end + # Next step shake + v_true = nextv() + v_cumul .+= v_true + end +end + +function imgen_square(sz) + return ImGen(curry(generate_square, sz), sz, 1, 1, "square$(sz[1])x$(sz[2])") +end + +################ +# Shake a photo +################ + +function generate_shake_image(im, sz, + :: Type{DisplacementConstant}, + datachannel :: Channel{OnlineData{DisplacementConstant}}, + params :: NamedTuple) + + # Restart the seed to enable comparison across predictors + Random.seed!(rng,9182737465) + + nextv = shake(params) + v_true = nextv() + v_cumul = copy(v_true) + + while true + # Extract subwindow of original image and add noise + b_true = zeros(sz...) + extract_subimage!(b_true, im, v_cumul; threads=true) + b = b_true .+ params.noise_level.*randn(rng,sz...) + v = v_true.*(1.0 .+ params.shake_noise_level.*randn(rng,size(v_true)...)) + # Pass data to iteration routine + data = OnlineData{DisplacementConstant}(b_true, b, v, v_true, v_cumul) + if !put_unless_closed!(datachannel, data) + return + end + # Next step shake + v_true = nextv() + v_cumul .+= v_true + end +end + +function imgen_shake(imname, sz) + im = Float64.(Gray.(TestImages.testimage(imname))) + dynrange = maximum(im) + return ImGen(curry(generate_shake_image, im, sz), sz, 1, dynrange, + "$(imname)$(sz[1])x$(sz[2])") +end + + + +######################################################################## +# PETscan +######################################################################## +function generate_radon(im, sz, + :: Type{DisplacementConstant}, + datachannel :: Channel{PetOnlineData{DisplacementConstant}}, + params :: NamedTuple) + + # Restart the seed to enable comparison across predictors + Random.seed!(rng,9182737465) + + nextv = shake(params) + v_true = nextv() + v_cumul = copy(v_true) + + S_true = generate_radonmask(params) + + theta_true = rotatebytheta(params) + theta_cumul = copy(theta_true) + + while true + # Define the transformation matrix + center_point = (sz[1]/2 + v_true[1], sz[2]/2 + v_true[2]) + tform = recenter(RotMatrix(theta_cumul), center_point) + + # Apply the transformation to the image using warp + b_true = copy(warp(im, tform, axes(im), fillvalue=Flat())) + + v = v_true.*(1.0 .+ params.shake_noise_level.*randn(rng,size(v_true)...)) + theta = theta_true*(1.0 + params.rotation_noise_level.*randn(rng)) + + # Generate the true and noisy sinograms + sinogram_true = zeros(params.radondims...) + sinogram_true .*= params.scale + radon!(sinogram_true, b_true) + sinogram_noisy = copy(sinogram_true) + + for i=1:params.radondims[1], j=1:params.radondims[2] + sinogram_noisy[i, j] += pois_rand(rng,params.noise_level) + end + + # Pass data to iteration routine + data = PetOnlineData{DisplacementConstant}(b_true, sinogram_true, sinogram_noisy, v, v_true, v_cumul, [theta, theta_cumul], S_true) + if !put_unless_closed!(datachannel, data) + return + end + + # Next step shake + v_true = nextv() + v_cumul .+= v_true + # Next theta + theta_true = rotatebytheta(params) + theta_cumul += theta_true + # Next sinogram mask + S_true = generate_radonmask(params) + end +end + + +function imgen_shepplogan_radon(origsize,sz) + im = convert(Array{Float64},TestImages.shepp_logan(origsize, highContrast=true)) + dynrange = maximum(im) + return ImGen(curry(generate_radon, im, sz), sz, 1, dynrange, "shepplogan$(sz[1])x$(sz[2])") +end + + +end # Module
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/src/PET/OpticalFlow.jl Fri Apr 19 17:00:37 2024 +0300 @@ -0,0 +1,431 @@ +################################ +# Code relevant to optical flow +################################ + +__precompile__() + +module OpticalFlow + +using AlgTools.Util +using ImageTools.Gradient +import ImageTools.Translate +using ImageTools.ImFilter + +# using ImageTransformations +# using Images, CoordinateTransformations, Rotations, OffsetArrays +# using Interpolations + +import Images: center, warp +import CoordinateTransformations: recenter +import Rotations: RotMatrix +import Interpolations: Flat + +########## +# Exports +########## + +export flow!, + pdflow!, + flow_grad!, + flow_interp!, + estimate_Λ², + estimate_linear_Λ², + pointwise_gradiprod_2d!, + pointwise_gradiprod_2dᵀ!, + horn_schunck_reg_prox!, + horn_schunck_reg_prox_op!, + mldivide_step_plus_sym2x2!, + linearised_optical_flow_error, + Image, AbstractImage, ImageSize, + Gradient, Displacement, + DisplacementFull, DisplacementConstant, + HornSchunckData, + filter_hs, + petpdflow! + +############################################### +# Types (several imported from ImageTools.Translate) +############################################### + +Image = Translate.Image +AbstractImage = AbstractArray{Float64,2} +Displacement = Translate.Displacement +DisplacementFull = Translate.DisplacementFull +DisplacementConstant = Translate.DisplacementConstant +Gradient = Array{Float64,3} +ImageSize = Tuple{Int64,Int64} + +################################# +# Displacement field based flow +################################# + +function flow_interp!(x::AbstractImage, u::Displacement, tmp::AbstractImage; + threads = false) + tmp .= x + Translate.translate_image!(x, tmp, u; threads=threads) +end + +function flow_interp!(x::AbstractImage, u::Displacement; + threads = false) + tmp = copy(x) + Translate.translate_image!(x, tmp, u; threads=threads) +end + +flow! = flow_interp! + +function pdflow!(x, Δx, y, Δy, u, dual_flow; threads=:none) + if dual_flow + #flow!((x, @view(y[1, :, :]), @view(y[2, :, :])), diffu, + # (Δx, @view(Δy[1, :, :]), @view(Δy[2, :, :]))) + @backgroundif (threads==:outer) begin + flow!(x, u, Δx; threads=(threads==:inner)) + end begin + flow!(@view(y[1, :, :]), u, @view(Δy[1, :, :]); threads=(threads==:inner)) + flow!(@view(y[2, :, :]), u, @view(Δy[2, :, :]); threads=(threads==:inner)) + end + else + flow!(x, u, Δx) + end +end + +function pdflow!(x, Δx, y, Δy, z, Δz, u, dual_flow; threads=:none) + if dual_flow + @backgroundif (threads==:outer) begin + flow!(x, u, Δx; threads=(threads==:inner)) + flow!(z, u, Δz; threads=(threads==:inner)) + end begin + flow!(@view(y[1, :, :]), u, @view(Δy[1, :, :]); threads=(threads==:inner)) + flow!(@view(y[2, :, :]), u, @view(Δy[2, :, :]); threads=(threads==:inner)) + end + else + flow!(x, u, Δx; threads=(threads==:inner)) + flow!(z, u, Δz; threads=(threads==:inner)) + end +end + +# Additional method for Greedy +function pdflow!(x, Δx, y, Δy, u; threads=:none) + @assert(size(u)==(2,)) + Δx .= x + Δy .= y + flow!(x, u; threads=(threads==:inner)) + Dxx = similar(Δy) + DΔx = similar(Δy) + ∇₂!(Dxx, x) + ∇₂!(DΔx, Δx) + inds = abs.(Dxx) .≤ 1e-1 + Dxx[inds] .= 1 + DΔx[inds] .= 1 + y .= y.* DΔx ./ Dxx +end + +# Additional method for Rotation +function pdflow!(x, Δx, y, u; threads=:none) + @assert(size(u)==(2,)) + Δx .= x + flow!(x, u; threads=(threads==:inner)) + + (m,n) = size(x) + dx = similar(y) + dx_banana = similar(y) + ∇₂!(dx, Δx) + ∇₂!(dx_banana, x) + + for i=1:m + for j=1:n + ndx = @views sum(dx[:, i, j].^2) + ndx_banana = @views sum(dx_banana[:, i, j].^2) + if ndx > 1e-4 && ndx_banana > 1e-4 + A = dx[:, i, j] + B = dx_banana[:, i, j] + theta = atan(B[1] * A[2] - B[2] * A[1], B[1] * A[1] + B[2] * A[2]) # Oriented angle from A to B + cos_theta = cos(theta) + sin_theta = sin(theta) + a = cos_theta * y[1, i, j] - sin_theta * y[2, i, j] + b = sin_theta * y[1, i, j] + cos_theta * y[2, i, j] + y[1, i, j] = a + y[2, i, j] = b + end + end + end +end + +# Additional method for Dual Scaling +function pdflow!(x, y, u; threads=:none) + @assert(size(u)==(2,)) + oldx = copy(x) + flow!(x, u; threads=(threads==:inner)) + C = similar(y) + cc = abs.(x-oldx) + cm = max(1e-12,maximum(cc)) + c = 1 .* (1 .- cc./ cm) .^(10) + C[1,:,:] .= c + C[2,:,:] .= c + y .= C.*y +end + + +########################## +# PET +########################## +function petflow_interp!(x::AbstractImage, tmp::AbstractImage, u::DisplacementConstant, theta_known::DisplacementConstant; + threads = false) + tmp .= x + center_point = center(x) .+ u + tform = recenter(RotMatrix(theta_known[1]), center_point) + tmp = warp(x, tform, axes(x), fillvalue=Flat()) + x .= tmp +end + +petflow! = petflow_interp! + +function petpdflow!(x, Δx, y, Δy, u, theta_known, dual_flow; threads=:none) + if dual_flow + @backgroundif (threads==:outer) begin + petflow!(x, Δx, u, theta_known; threads=(threads==:inner)) + end begin + petflow!(@view(y[1, :, :]), @view(Δy[1, :, :]), u, theta_known; threads=(threads==:inner)) + petflow!(@view(y[2, :, :]), @view(Δy[2, :, :]), u, theta_known; threads=(threads==:inner)) + end + else + petflow!(x, Δx, u, theta_known) + end +end + +# Method for greedy predictor +function petpdflow!(x, Δx, y, Δy, u, theta_known, dual_flow, β; threads=:none) + oldx = copy(x) + center_point = center(x) .+ u + tform = recenter(RotMatrix(theta_known[1]), center_point) + Δx = warp(x, tform, axes(x), fillvalue=Flat()) + @. x = Δx + @. Δy = y + if dual_flow + Dxx = copy(Δy) + DΔx = copy(Δy) + ∇₂!(Dxx, x) + ∇₂!(DΔx, oldx) + inds = abs.(Dxx) .≤ β + Dxx[inds] .= 1 + DΔx[inds] .= 1 + y .= y.* DΔx ./ Dxx + end +end + +# Method for affine predictor +function petpdflow!(x, Δx, y, u, theta_known, dual_flow; threads=:none) + oldx = copy(x) + center_point = center(x) .+ u + tform = recenter(RotMatrix(theta_known[1]), center_point) + Δx = warp(x, tform, axes(x), fillvalue=Flat()) + @. x = Δx + C = similar(y) + cc = abs.(x-oldx) + if dual_flow + cm = max(1e-12,maximum(cc)) + c = 1 .* (1 .- cc./ cm) .^(10) + C[1,:,:] .= c + C[2,:,:] .= c + y .= C.*y + end +end + +# Method for rotation prediction (exploiting property of inverse rotation) +function petpdflow!(x, Δx, y, Δy, u, theta_known, dual_flow, β₁, β₂; threads=:none) + if dual_flow + @backgroundif (threads==:outer) begin + petflow!(x, Δx, u, theta_known; threads=(threads==:inner)) + end begin + petflow!(@view(y[1, :, :]), @view(Δy[1, :, :]), u, -theta_known; threads=(threads==:inner)) + petflow!(@view(y[2, :, :]), @view(Δy[2, :, :]), u, -theta_known; threads=(threads==:inner)) + end + else + petflow!(x, Δx, u, theta_known) + end +end + +########################## +# Linearised optical flow +########################## + +# ⟨⟨u, ∇b⟩⟩ +function pointwise_gradiprod_2d!(y::Image, vtmp::Gradient, + u::DisplacementFull, b::Image; + add = false) + ∇₂c!(vtmp, b) + + u′=reshape(u, (size(u, 1), prod(size(u)[2:end]))) + vtmp′=reshape(vtmp, (size(vtmp, 1), prod(size(vtmp)[2:end]))) + y′=reshape(y, prod(size(y))) + + if add + @simd for i = 1:length(y′) + @inbounds y′[i] += dot(@view(u′[:, i]), @view(vtmp′[:, i])) + end + else + @simd for i = 1:length(y′) + @inbounds y′[i] = dot(@view(u′[:, i]), @view(vtmp′[:, i])) + end + end +end + +function pointwise_gradiprod_2d!(y::Image, vtmp::Gradient, + u::DisplacementConstant, b::Image; + add = false) + ∇₂c!(vtmp, b) + + vtmp′=reshape(vtmp, (size(vtmp, 1), prod(size(vtmp)[2:end]))) + y′=reshape(y, prod(size(y))) + + if add + @simd for i = 1:length(y′) + @inbounds y′[i] += dot(u, @view(vtmp′[:, i])) + end + else + @simd for i = 1:length(y′) + @inbounds y′[i] = dot(u, @view(vtmp′[:, i])) + end + end +end + +# ∇b ⋅ y +function pointwise_gradiprod_2dᵀ!(u::DisplacementFull, y::Image, b::Image) + ∇₂c!(u, b) + + u′=reshape(u, (size(u, 1), prod(size(u)[2:end]))) + y′=reshape(y, prod(size(y))) + + @simd for i=1:length(y′) + @inbounds @. u′[:, i] *= y′[i] + end +end + +function pointwise_gradiprod_2dᵀ!(u::DisplacementConstant, y::Image, b::Image) + @assert(size(y)==size(b) && size(u)==(2,)) + u .= 0 + ∇₂cfold!(b, nothing) do g, st, (i, j) + @inbounds u .+= g.*y[i, j] + return st + end + # Reweight to be with respect to 𝟙^*𝟙 inner product. + u ./= prod(size(b)) +end + +mutable struct ConstantDisplacementHornSchunckData + M₀::Array{Float64,2} + z::Array{Float64,1} + Mv::Array{Float64,2} + av::Array{Float64,1} + cv::Float64 + + function ConstantDisplacementHornSchunckData() + return new(zeros(2, 2), zeros(2), zeros(2,2), zeros(2), 0) + end +end + +# For DisplacementConstant, for the simple prox step +# +# (1) argmin_u 1/(2τ)|u-ũ|^2 + (θ/2)|b⁺-b+<<u-ŭ,∇b>>|^2 + (λ/2)|u-ŭ|^2, +# +# construct matrix M₀ and vector z such that we can solve u from +# +# (2) (I/τ+M₀)u = M₀ŭ + ũ/τ - z +# +# Note that the problem +# +# argmin_u 1/(2τ)|u-ũ|^2 + (θ/2)|b⁺-b+<<u-ŭ,∇b>>|^2 + (λ/2)|u-ŭ|^2 +# + (θ/2)|b⁺⁺-b⁺+<<uʹ-u,∇b⁺>>|^2 + (λ/2)|u-uʹ|^2 +# +# has with respect to u the system +# +# (I/τ+M₀+M₀ʹ)u = M₀ŭ + M₀ʹuʹ + ũ/τ - z + zʹ, +# +# where the primed variables correspond to (2) for (1) for uʹ in place of u: +# +# argmin_uʹ 1/(2τ)|uʹ-ũʹ|^2 + (θ/2)|b⁺⁺-b⁺+<<uʹ-u,∇b⁺>>|^2 + (λ/2)|uʹ-u|^2 +# +function horn_schunck_reg_prox_op!(hs::ConstantDisplacementHornSchunckData, + bnext::Image, b::Image, θ, λ, T) + @assert(size(b)==size(bnext)) + w = prod(size(b)) + z = hs.z + cv = 0 + # Factors of symmetric matrix [a c; c d] + a, c, d = 0.0, 0.0, 0.0 + # This used to use ∇₂cfold but it is faster to allocate temporary + # storage for the full gradient due to probably better memory and SIMD + # instruction usage. + g = zeros(2, size(b)...) + ∇₂c!(g, b) + @inbounds for i=1:size(b, 1) + for j=1:size(b, 2) + δ = bnext[i,j]-b[i,j] + @. z += g[:,i,j]*δ + cv += δ*δ + a += g[1,i,j]*g[1,i,j] + c += g[1,i,j]*g[2,i,j] + d += g[2,i,j]*g[2,i,j] + end + end + w₀ = λ + w₂ = θ/w + aʹ = w₀ + w₂*a + cʹ = w₂*c + dʹ = w₀ + w₂*d + hs.M₀ .= [aʹ cʹ; cʹ dʹ] + hs.Mv .= [w*λ+θ*a θ*c; θ*c w*λ+θ*d] + hs.cv = cv*θ + hs.av .= hs.z.*θ + hs.z .*= w₂/T +end + +# Solve the 2D system (I/τ+M₀)u = z +@inline function mldivide_step_plus_sym2x2!(u, M₀, z, τ) + a = 1/τ+M₀[1, 1] + c = M₀[1, 2] + d = 1/τ+M₀[2, 2] + u .= ([d -c; -c a]*z)./(a*d-c*c) +end + +function horn_schunck_reg_prox!(u::DisplacementConstant, bnext::Image, b::Image, + θ, λ, T, τ) + hs=ConstantDisplacementHornSchunckData() + horn_schunck_reg_prox_op!(hs, bnext, b, θ, λ, T) + mldivide_step_plus_sym2x2!(u, hs.M₀, (u./τ)-hs.z, τ) +end + +function flow_grad!(x::Image, vtmp::Gradient, u::Displacement; δ=nothing) + if !isnothing(δ) + u = δ.*u + end + pointwise_gradiprod_2d!(x, vtmp, u, x; add=true) +end + +# Error b-b_prev+⟨⟨u, ∇b⟩⟩ for Horn–Schunck type penalisation +function linearised_optical_flow_error(u::Displacement, b::Image, b_prev::Image) + imdim = size(b) + vtmp = zeros(2, imdim...) + tmp = b-b_prev + pointwise_gradiprod_2d!(tmp, vtmp, u, b_prev; add=true) + return tmp +end + +############################################## +# Helper to smooth data for Horn–Schunck term +############################################## + +function filter_hs(b, b_next, b_next_filt, kernel) + if kernel==nothing + f = x -> x + else + f = x -> simple_imfilter(x, kernel; threads=true) + end + + # We already filtered b in the previous step (b_next in that step) + b_filt = b_next_filt==nothing ? f(b) : b_next_filt + b_next_filt = f(b_next) + + return b_filt, b_next_filt +end + +end # Module
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/src/PET/PET.jl Fri Apr 19 17:00:37 2024 +0300 @@ -0,0 +1,506 @@ +__precompile__() + +module PET + +######################## +# Load external modules +######################## + +using Printf +using FileIO +#using JLD2 +using Setfield +using ImageQualityIndexes: assess_psnr, assess_ssim +using DelimitedFiles +import GR + +using AlgTools.Util +using AlgTools.StructTools +using AlgTools.LinkedLists +using AlgTools.Comms +using ImageTools.Visualise: secs_ns, grayimg, do_visualise +using ImageTools.ImFilter: gaussian + +# For PET +#using Colors +#using ColorTypes +using ColorSchemes +#using PerceptualColourMaps + +##################### +# Load local modules +##################### +include("OpticalFlow.jl") +include("Radon.jl") +include("ImGenerate.jl") +include("AlgorithmDualScaling.jl") +include("AlgorithmGreedy.jl") +include("AlgorithmNoPrediction.jl") +include("AlgorithmPrimalOnly.jl") +include("AlgorithmProximal.jl") +include("AlgorithmRotation.jl") + +using .Radon +using .ImGenerate +using .OpticalFlow +using .AlgorithmDualScaling +using .AlgorithmGreedy +using .AlgorithmNoPrediction +using .AlgorithmPrimalOnly +using .AlgorithmProximal +using .AlgorithmRotation + + +############## +# Our exports +############## + +export demo_pet1, demo_pet2, demo_pet3, demo_pet4, demo_pet5, demo_pet6, batchrun_pet + +################################### +# Parameterisation and experiments +################################### + +struct Experiment + mod :: Module + DisplacementT :: Type + imgen :: ImGen + params :: NamedTuple +end + +function Base.show(io::IO, e::Experiment) + displacementname(::Type{DisplacementFull}) = "DisplacementFull" + displacementname(::Type{DisplacementConstant}) = "DisplacementConstant" + print(io, " + mod: $(e.mod) + DisplacementT: $(displacementname(e.DisplacementT)) + imgen: $(e.imgen.name) $(e.imgen.dim[1])×$(e.imgen.dim[2]) + params: $(e.params ⬿ (kernel = "(not shown)",)) + ") +end + +const default_save_prefix="img/" + +const default_params = ( + ρ = 0, + verbose_iter = 100, + maxiter = 4000, + save_results = true, + save_images = true, + save_images_iters = Set([100, 200, 300, 500, 800, 1000, + 1300, 1500, 1800, 2000, 4000]), + pixelwise_displacement=false, + dual_flow = true, + prox_predict = true, + handle_interrupt = true, + init = :zero, + plot_movement = false, +) + +const p_known₀_pet = ( + noise_level = 0.5, + shake_noise_level = 0.1, + shake = 1.0, + rotation_factor = 0.15, + rotation_noise_level = 0.0025, + α = 1.0, + ρ̃₀ = 1.0, + σ̃₀ = 0.0001, + δ = 0.9, + σ₀ = 0.5, + τ₀ = 0.9, + λ = 1, + origsize = 256, + radondims = [128, 64], + sz = (256,256), + scale = 1, + c = 1.0, + sino_sparsity = 0.5, + L = 300.0, + L_experiment = false, +) + +const petscan = imgen_shepplogan_radon(p_known₀_pet.origsize, p_known₀_pet.sz) + +const pet_experiments_pdps_known = ( + Experiment(AlgorithmDualScaling, DisplacementConstant, petscan, + p_known₀_pet), + Experiment(AlgorithmGreedy, DisplacementConstant, petscan, + p_known₀_pet), + Experiment(AlgorithmNoPrediction, DisplacementConstant, petscan, + p_known₀_pet), + Experiment(AlgorithmPrimalOnly, DisplacementConstant, petscan, + p_known₀_pet), + Experiment(AlgorithmProximal, DisplacementConstant, petscan, + p_known₀_pet ⬿ (phantom_ρ = 100,)), + Experiment(AlgorithmRotation, DisplacementConstant, petscan, + p_known₀_pet), +) + +const pet_experiments_all = Iterators.flatten(( + pet_experiments_pdps_known, +)) + +################ +# Log +################ + +struct LogEntry <: IterableStruct + iter :: Int + time :: Float64 + function_value :: Float64 + #v_cumul_true_y :: Float64 + #v_cumul_true_x :: Float64 + #v_cumul_est_y :: Float64 + #v_cumul_est_x :: Float64 + psnr :: Float64 + ssim :: Float64 + #psnr_data :: Float64 + #ssim_data :: Float64 +end + +struct LogEntryHiFi <: IterableStruct + iter :: Int + v_cumul_true_y :: Float64 + v_cumul_true_x :: Float64 +end + +############### +# Main routine +############### + +struct State + vis :: Union{Channel,Bool,Nothing} + start_time :: Union{Real,Nothing} + wasted_time :: Real + log :: LinkedList{LogEntry} + log_hifi :: LinkedList{LogEntryHiFi} + aborted :: Bool +end + +function name(e::Experiment, p) + ig = e.imgen + # return "$(ig.name)_$(e.mod.identifier)_$(@sprintf "%x" hash(p))" + return "$(ig.name)_$(e.mod.identifier)_$(Int64(10000*p.σ₀))_$(Int64(10000*p.τ₀))" +end + +function write_tex(texfile, e_params) + open(texfile, "w") do io + wp = (n, v) -> println(io, "\\def\\EXPPARAM$(n){$(v)}") + wf = (n, s) -> if isdefined(e_params, s) + wp(n, getfield(e_params, s)) + end + wf("alpha", :α) + wf("sigmazero", :σ₀) + wf("tauzero", :τ₀) + wf("tildetauzero", :τ̃₀) + wf("delta", :δ) + wf("lambda", :λ) + wf("theta", :θ) + wf("maxiter", :maxiter) + wf("noiselevel", :noise_level) + wf("shakenoiselevel", :shake_noise_level) + wf("shake", :shake) + wf("timestep", :timestep) + wf("displacementcount", :displacementcount) + wf("phantomrho", :phantom_ρ) + if isdefined(e_params, :σ₀) + wp("sigma", (e_params.σ₀ == 1 ? "" : "$(e_params.σ₀)") * "\\sigma_{\\max}") + end + end +end + +function run_experiments(;visualise=true, + recalculate=true, + experiments, + save_prefix=default_save_prefix, + fullscreen=false, + kwargs...) + + # Create visualisation + if visualise + rc = Channel(1) + visproc = Threads.@spawn bg_visualise_enhanced(rc, fullscreen=fullscreen) + bind(rc, visproc) + vis = rc + else + vis = false + end + + # Run all experiments + for e ∈ experiments + + # Parameters for this experiment + e_params = default_params ⬿ e.params ⬿ kwargs + ename = name(e, e_params) + e_params = e_params ⬿ (save_prefix = save_prefix * ename, + dynrange = e.imgen.dynrange, + Λ = e.imgen.Λ) + + if recalculate || !isfile(e_params.save_prefix * ".txt") + println("Running experiment \"$(ename)\"") + + # Start data generation task + datachannel = Channel{PetOnlineData{e.DisplacementT}}(2) + gentask = Threads.@spawn e.imgen.f(e.DisplacementT, datachannel, e_params) + bind(datachannel, gentask) + + # Run algorithm + iterate = curry(iterate_visualise, datachannel, + State(vis, nothing, 0.0, nothing, nothing, false)) + + x, y, st = e.mod.solve(e.DisplacementT; + dim=e.imgen.dim, + iterate=iterate, + params=e_params) + + # Clear non-saveable things + st = @set st.vis = nothing + + println("Wasted_time: $(st.wasted_time)s") + + if e_params.save_results + println("Saving " * e_params.save_prefix * "(.txt,_hifi.txt,_params.tex)") + + perffile = e_params.save_prefix * ".txt" + hififile = e_params.save_prefix * "_hifi.txt" + texfile = e_params.save_prefix * "_params.tex" + # datafile = e_params.save_prefix * ".jld2" + + write_log(perffile, st.log, "# params = $(e_params)\n") + write_log(hififile, st.log_hifi, "# params = $(e_params)\n") + write_tex(texfile, e_params) + # @save datafile x y st params + end + + close(datachannel) + wait(gentask) + + if st.aborted + break + end + else + println("Skipping already computed experiment \"$(ename)\"") + # texfile = e_params.save_prefix * "_params.tex" + # write_tex(texfile, e_params) + end + end + + if visualise + # Tell subprocess to finish, and wait + put!(rc, nothing) + close(rc) + wait(visproc) + end + + return +end + +####################### +# Demos and batch runs +####################### + +function demo(experiment; kwargs...) + run_experiments(;experiments=(experiment,), + save_results=false, + save_images=false, + visualise=true, + recalculate=true, + verbose_iter=1, + fullscreen=true, + kwargs...) +end + +demo_pet1 = () -> demo(pet_experiments_pdps_known[1]) # Dual scaling +demo_pet2 = () -> demo(pet_experiments_pdps_known[2]) # Greedy +demo_pet3 = () -> demo(pet_experiments_pdps_known[3]) # No Prediction +demo_pet4 = () -> demo(pet_experiments_pdps_known[4]) # Primal only +demo_pet5 = () -> demo(pet_experiments_pdps_known[5]) # Proximal (old) +demo_pet6 = () -> demo(pet_experiments_pdps_known[6]) # Rotation + + +function batchrun_pet(kwargs...) + run_experiments(;experiments=pet_experiments_all, + save_results=true, + save_images=true, + visualise=false, + recalculate=false, + kwargs...) +end + +###################################################### +# Iterator that does visualisation and log collection +###################################################### + +function rescale(arr, new_range) + old_min = minimum(arr) + old_max = maximum(arr) + + scale_factor = (new_range[2] - new_range[1]) / (old_max - old_min) + + scaled_arr = new_range[1] .+ (arr .- old_min) * scale_factor + + return scaled_arr +end + +function iterate_visualise(datachannel::Channel{PetOnlineData{DisplacementT}}, + st :: State, + step :: Function, + params :: NamedTuple) where DisplacementT + try + sc = nothing + + d = take!(datachannel) + + for iter=1:params.maxiter + dnext = take!(datachannel) + st = step(d.sinogram_noisy, d.v, d.theta, d.b_true, d.S) do calc_objective + stn = st + + if isnothing(stn.start_time) + # The Julia precompiler is a miserable joke, apparently not crossing module + # boundaries, so only start timing after the first iteration. + stn = @set stn.start_time=secs_ns() + end + + verb = params.verbose_iter!=0 && mod(iter, params.verbose_iter) == 0 + + # Normalise movement to image dimensions so + # our TikZ plotting code doesn't need to know + # the image pixel size. + sc = 1.0./maximum(size(d.b_true)) + + if verb || iter ≤ 20 || (iter ≤ 200 && mod(iter, 10) == 0) + verb_start = secs_ns() + tm = verb_start - stn.start_time - stn.wasted_time + value, x, v, vhist = calc_objective() + + entry = LogEntry(iter, tm, value, + #sc*d.v_cumul_true[1], + #sc*d.v_cumul_true[2], + #sc*v[1], sc*v[2], + assess_psnr(x, d.b_true), + assess_ssim(x, d.b_true), + #assess_psnr(d.b_noisy, d.b_true), + #assess_ssim(d.b_noisy, d.b_true) + ) + + # (**) Collect a singly-linked list of log to avoid array resizing + # while iterating + stn = @set stn.log=LinkedListEntry(entry, stn.log) + + if !isnothing(vhist) + vhist=vhist.*sc + end + + if verb + @printf("%d/%d J=%f, PSNR=%f, SSIM=%f, avg. FPS=%f\n", + iter, params.maxiter, value, entry.psnr, + entry.ssim, entry.iter/entry.time) + if isa(stn.vis, Channel) + put_onlylatest!(stn.vis, ((rescale(backproject!(d.b_true,d.sinogram_noisy),(0.0,params.dynrange)), x), + params.plot_movement, + stn.log, vhist)) + + end + end + + if params.save_images && (!haskey(params, :save_images_iters) + || iter ∈ params.save_images_iters) + fn = (t, ext) -> "$(params.save_prefix)_$(t)_frame$(iter).$(ext)" + normalise = (data) -> data./maximum(data) + # save(File(format"PNG", fn("true", "png")), mapped_img(d.b_true, ColorSchemes.cmyk.colors[1:end])) + # save(File(format"PNG", fn("true_sinogram", "png")), mapped_img(normalise(d.sinogram_true), ColorSchemes.cmyk.colors[1:end])) + # save(File(format"PNG", fn("data_sinogram", "png")), mapped_img(normalise(d.S.*d.sinogram_noisy), ColorSchemes.cmyk.colors[1:end])) + save(File(format"PNG", fn("reco", "png")), mapped_img(x, ColorSchemes.cmyk.colors[1:end])) + if !isnothing(vhist) + open(fn("movement", "txt"), "w") do io + writedlm(io, ["est_y" "est_x"]) + writedlm(io, vhist) + end + end + end + + stn = @set stn.wasted_time += (secs_ns() - verb_start) + + return stn + end + + hifientry = LogEntryHiFi(iter, sc*d.v_cumul_true[1], sc*d.v_cumul_true[2]) + st = @set st.log_hifi=LinkedListEntry(hifientry, st.log_hifi) + + return st + end + d=dnext + end + catch ex + if params.handle_interrupt && isa(ex, InterruptException) + # If SIGINT is received (user pressed ^C), terminate computations, + # returning current status. Effectively, we do not call `step()` again, + # ending the iterations, but letting the algorithm finish up. + # Assuming (**) above occurs atomically, `st.log` should be valid, but + # any results returned by the algorithm itself may be partial, as for + # reasons of efficiency we do *not* store results of an iteration until + # the next iteration is finished. + printstyled("\rUser interrupt—finishing up.\n", bold=true, color=202) + st = @set st.aborted = true + else + rethrow(ex) + end + end + + return st +end + +function bg_visualise_enhanced(rc; fullscreen=false) + process_channel(rc) do d + imgs, plot_movement, log, vhist = d + do_visualise(imgs, refresh=false, fullscreen=fullscreen) + # Overlay movement + GR.settextcolorind(5) + GR.setcharheight(0.015) + GR.settextpath(GR.TEXT_PATH_RIGHT) + tx, ty = GR.wctondc(0, 1) + GR.text(tx, ty, @sprintf "FPS %.1f, SSIM %.2f, PSNR %.1f" (log.value.iter/log.value.time) log.value.ssim log.value.psnr) + if plot_movement + sc=1.0 + p=unfold_linked_list(log) + x=map(e -> 1.5+sc*e.v_cumul_true_x, p) + y=map(e -> 0.5+sc*e.v_cumul_true_y, p) + GR.setlinewidth(2) + GR.setlinecolorind(2) + GR.polyline(x, y) + x=map(e -> 1.5+sc*e.v_cumul_est_x, p) + y=map(e -> 0.5+sc*e.v_cumul_est_y, p) + GR.setlinecolorind(3) + GR.polyline(x, y) + if vhist != nothing + GR.setlinecolorind(4) + x=map(v -> 1.5+sc*v, vhist[:,2]) + y=map(v -> 0.5+sc*v, vhist[:,1]) + GR.polyline(x, y) + end + end + GR.updatews() + end +end + +# Clip image values to allowed range +clip = x -> min(max(x, 0.0), 1.0) + +# Apply a colourmap (vector of RGB objects) to raw image data +function mapped_img(im, cmap) + l = length(cmap) + apply = t -> cmap[1+round(UInt16, clip(t) * (l-1))] + return apply.(im) +end + +############### +# Precompiling +############### + +# precompile(Tuple{typeof(GR.drawimage), Float64, Float64, Float64, Float64, Int64, Int64, Array{UInt32, 2}}) +# precompile(Tuple{Type{Plots.Plot{T} where T<:RecipesBase.AbstractBackend}, Plots.GRBackend, Int64, Base.Dict{Symbol, Any}, Base.Dict{Symbol, Any}, Array{Plots.Series, 1}, Nothing, Array{Plots.Subplot{T} where T<:RecipesBase.AbstractBackend, 1}, Base.Dict{Any, Plots.Subplot{T} where T<:RecipesBase.AbstractBackend}, Plots.EmptyLayout, Array{Plots.Subplot{T} where T<:RecipesBase.AbstractBackend, 1}, Bool}) +# precompile(Tuple{typeof(Plots._plot!), Plots.Plot{Plots.GRBackend}, Base.Dict{Symbol, Any}, Tuple{Array{ColorTypes.Gray{Float64}, 2}}}) + +end # Module
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/src/PET/Radon.jl Fri Apr 19 17:00:37 2024 +0300 @@ -0,0 +1,85 @@ +# From https://github.com/LMescheder/Radon.jl/blob/master/src/Radon.jl +# Updated to modern Julia. Fixed radon! to zero out `radonim`. + +__precompile__() + + +module Radon +export radon!, backproject! +# package code goes here + +function radon(im::AbstractArray{T, 2}) where T<:AbstractFloat + Nx, Ny = size(im) + radonim = similar(im, min(Nx, Ny), max(Nx, Ny)) + radon!(radonim, im) + radonim +end + +function radon!(radonim::AbstractArray{T, 2}, im::AbstractArray{T, 2}) where T<:AbstractFloat + Nφ, Ns = size(radonim) + Nx, Ny = size(im) + Ldiag = hypot(Nx, Ny) + K = round(Int, Ldiag) + + dLz = Ldiag/(K-1) + dLs = Ldiag/(Ns-1) + dφ = π/Nφ + + @inbounds for is = 1:Ns, iφ = 1:Nφ + φ = (iφ-1)*dφ + s = -Ldiag/2 + (is-1)*dLs + tx = cos(φ) + ty = sin(φ) + + radonim[iφ, is] = zero(T) + for k=1:K + z = -Ldiag/2 + (k-1)*dLz + x = round(Int, Nx/2 + z*ty + s*tx) + y = round(Int, Ny/2 - z*tx + s*ty) + + if 1 <= x <= Nx && 1 <= y <= Ny + radonim[iφ, is] += im[x, y]*dLz + end + end + end + radonim +end + +function backproject(radonim::AbstractArray{T, 2}) where T<:AbstractFloat + Nφ, Ns = size(radonim) + im = similar(radonim, Ns, Ns) + backproject!(im, radonim) + im +end + +function backproject!(im::AbstractArray{T, 2}, radonim::AbstractArray{T, 2}) where T<:AbstractFloat + Nφ, Ns = size(radonim) + Nx, Ny = size(im) + Ldiag = hypot(Nx, Ny) + K = round(Int, Ldiag) + + dLz = Ldiag/(K-1) + dLs = Ldiag/(Ns-1) + dφ = π/Nφ + + im .= 0. + @inbounds for is = 1:Ns, iφ = 1:Nφ + φ = (iφ-1)*dφ + s = -Ldiag/2 + (is-1)*dLs + tx = cos(φ) + ty = sin(φ) + + for k=1:K + z = -Ldiag/2 + (k-1)*dLz + x = round(Int, Nx/2 + z*ty + s*tx) + y = round(Int, Ny/2 - z*tx + s*ty) + + if 1 <= x <= Nx && 1 <= y <= Ny + im[x, y] += radonim[iφ, is]*dLz + end + end + end + im +end + +end # module
--- a/src/PredictPDPS.jl Thu Apr 18 11:31:32 2024 +0300 +++ b/src/PredictPDPS.jl Fri Apr 19 17:00:37 2024 +0300 @@ -30,6 +30,7 @@ ##################### include("OpticalFlow.jl") +include("Radon.jl") include("ImGenerate.jl") include("Algorithm.jl") include("AlgorithmBoth.jl") @@ -40,6 +41,7 @@ include("AlgorithmFB.jl") include("AlgorithmFBDual.jl") + # Additional include("AlgorithmProximal.jl") include("AlgorithmGreedy.jl") @@ -47,6 +49,7 @@ include("AlgorithmNoPrediction.jl") include("AlgorithmPrimalOnly.jl") include("AlgorithmDualScaling.jl") +include("PET/PET.jl") import .Algorithm, .AlgorithmBoth, @@ -61,10 +64,14 @@ .AlgorithmRotation, .AlgorithmNoPrediction, .AlgorithmPrimalOnly, - .AlgorithmDualScaling + .AlgorithmDualScaling, + .PET using .ImGenerate using .OpticalFlow: DisplacementFull, DisplacementConstant +using .PET: demo_pet1, demo_pet2, demo_pet3, + demo_pet4, demo_pet5, demo_pet6, + batchrun_pet ############## # Our exports @@ -72,23 +79,15 @@ export run_experiments, batchrun_article, - demo_known1, - demo_known2, - demo_known3, - demo_known4, - demo_known5, - demo_known6, - demo_unknown1, - demo_unknown2, - demo_unknown3, - demo_predictor1, - demo_predictor2, - demo_predictor3, - demo_predictor4, - demo_predictor5, - demo_predictor6, - batchrun_predictors - + demo_known1, demo_known2, demo_known3, + demo_unknown1,demo_unknown2,demo_unknown3, + batchrun_denoising, + demo_denoising1, demo_denoising2, demo_denoising3, + demo_denoising4, demo_denoising5, demo_denoising6, + demo_pet1, demo_pet2, demo_pet3, + demo_pet4, demo_pet5, demo_pet6, + batchrun_pet, batchrun_predictors + ################################### # Parameterisation and experiments ################################### @@ -137,13 +136,13 @@ const p_known₀ = ( noise_level = 0.5, - shake_noise_level = 0.025, - shake = 3.0, - α = 1.0, - ρ̃₀ = 1.0, - σ̃₀ = 1.0, + shake_noise_level = 0.05, + shake = 2, + α = 1, + ρ̃₀ = 1, + σ̃₀ = 1, δ = 0.9, - σ₀ = 0.1, + σ₀ = 1, τ₀ = 0.01, ) @@ -164,9 +163,21 @@ τ₀ = 0.01, ) +const p_known₀_denoising = ( + noise_level = 0.5, + shake_noise_level = 0.025, + shake = 3.0, + α = 1.0, + ρ̃₀ = 1.0, + σ̃₀ = 1.0, + δ = 0.9, + σ₀ = 0.1, + τ₀ = 0.01, +) + const experiments_pdps_known = ( Experiment(Algorithm, DisplacementConstant, lighthouse, - p_known₀ ⬿ (phantom_ρ = 0,)), + p_known₀_denoising ⬿ (phantom_ρ = 0,)), Experiment(Algorithm, DisplacementConstant, lighthouse, p_known₀ ⬿ (phantom_ρ = 100,)), Experiment(Algorithm, DisplacementConstant, square, @@ -187,32 +198,30 @@ p_known₀ ⬿ (τ̃₀=0.9, fb_inner_iterations = 10)), ) -const predictor_experiments_pdps_known = ( - Experiment(AlgorithmProximal, DisplacementConstant, lighthouse, - p_known₀ ⬿ (phantom_ρ = 100,)), - Experiment(AlgorithmRotation, DisplacementConstant, lighthouse, - p_known₀), - Experiment(AlgorithmGreedy, DisplacementConstant, lighthouse, - p_known₀), - Experiment(AlgorithmPrimalOnly, DisplacementConstant, lighthouse, - p_known₀), - Experiment(AlgorithmNoPrediction, DisplacementConstant, lighthouse, - p_known₀), - Experiment(AlgorithmDualScaling, DisplacementConstant, lighthouse, - p_known₀), -) - -const predictor_experiments_all = Iterators.flatten(( - predictor_experiments_pdps_known, -)) - const experiments_all = Iterators.flatten(( experiments_pdps_known, experiments_pdps_unknown_multi, experiments_fb_known )) +const denoising_experiments_pdps_known = ( + Experiment(AlgorithmDualScaling, DisplacementConstant, lighthouse, + p_known₀_denoising), + Experiment(AlgorithmGreedy, DisplacementConstant, lighthouse, + p_known₀_denoising), + Experiment(AlgorithmNoPrediction, DisplacementConstant, lighthouse, + p_known₀_denoising), + Experiment(AlgorithmPrimalOnly, DisplacementConstant, lighthouse, + p_known₀_denoising), + Experiment(AlgorithmProximal, DisplacementConstant, lighthouse, + p_known₀_denoising ⬿ (phantom_ρ = 100,)), + Experiment(AlgorithmRotation, DisplacementConstant, lighthouse, + p_known₀_denoising), +) +const denoising_experiments_all = Iterators.flatten(( + denoising_experiments_pdps_known, +)) ################ # Log @@ -388,17 +397,17 @@ demo_known2 = () -> demo(experiments_pdps_known[1]) demo_known3 = () -> demo(experiments_pdps_known[2]) -demo_predictor1 = () -> demo(predictor_experiments_pdps_known[1]) # Proximal -demo_predictor2 = () -> demo(predictor_experiments_pdps_known[2]) # Rotation -demo_predictor3 = () -> demo(predictor_experiments_pdps_known[3]) # Greedy -demo_predictor4 = () -> demo(predictor_experiments_pdps_known[4]) # PrimalOnly -demo_predictor5 = () -> demo(predictor_experiments_pdps_known[5]) # NoPrediction -demo_predictor6 = () -> demo(predictor_experiments_pdps_known[6]) # DualScaling - demo_unknown1 = () -> demo(experiments_pdps_unknown_multi[3], plot_movement=true) demo_unknown2 = () -> demo(experiments_pdps_unknown_multi[1], plot_movement=true) demo_unknown3 = () -> demo(experiments_pdps_unknown_multi[2], plot_movement=true) +demo_denoising1 = () -> demo(denoising_experiments_pdps_known[1]) # Dual scaling +demo_denoising2 = () -> demo(denoising_experiments_pdps_known[2]) # Greedy +demo_denoising3 = () -> demo(denoising_experiments_pdps_known[3]) # No Prediction +demo_denoising4 = () -> demo(denoising_experiments_pdps_known[4]) # Primal Only +demo_denoising5 = () -> demo(denoising_experiments_pdps_known[5]) # Proximal (old) +demo_denoising6 = () -> demo(denoising_experiments_pdps_known[6]) # Rotation + function batchrun_article(kwargs...) run_experiments(;experiments=experiments_all, save_results=true, @@ -408,8 +417,8 @@ kwargs...) end -function batchrun_predictors(kwargs...) - run_experiments(;experiments=predictor_experiments_all, +function batchrun_denoising(kwargs...) + run_experiments(;experiments=denoising_experiments_all, save_results=true, save_images=true, visualise=false, @@ -417,6 +426,12 @@ kwargs...) end + +function batchrun_predictors(kwargs...) + batchrun_denoising() + batchrun_pet() +end + ###################################################### # Iterator that does visualisation and log collection ######################################################
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/src/Radon.jl Fri Apr 19 17:00:37 2024 +0300 @@ -0,0 +1,85 @@ +# From https://github.com/LMescheder/Radon.jl/blob/master/src/Radon.jl +# Updated to modern Julia. Fixed radon! to zero out `radonim`. + +__precompile__() + + +module Radon +export radon!, backproject! +# package code goes here + +function radon(im::AbstractArray{T, 2}) where T<:AbstractFloat + Nx, Ny = size(im) + radonim = similar(im, min(Nx, Ny), max(Nx, Ny)) + radon!(radonim, im) + radonim +end + +function radon!(radonim::AbstractArray{T, 2}, im::AbstractArray{T, 2}) where T<:AbstractFloat + Nφ, Ns = size(radonim) + Nx, Ny = size(im) + Ldiag = hypot(Nx, Ny) + K = round(Int, Ldiag) + + dLz = Ldiag/(K-1) + dLs = Ldiag/(Ns-1) + dφ = π/Nφ + + @inbounds for is = 1:Ns, iφ = 1:Nφ + φ = (iφ-1)*dφ + s = -Ldiag/2 + (is-1)*dLs + tx = cos(φ) + ty = sin(φ) + + radonim[iφ, is] = zero(T) + for k=1:K + z = -Ldiag/2 + (k-1)*dLz + x = round(Int, Nx/2 + z*ty + s*tx) + y = round(Int, Ny/2 - z*tx + s*ty) + + if 1 <= x <= Nx && 1 <= y <= Ny + radonim[iφ, is] += im[x, y]*dLz + end + end + end + radonim +end + +function backproject(radonim::AbstractArray{T, 2}) where T<:AbstractFloat + Nφ, Ns = size(radonim) + im = similar(radonim, Ns, Ns) + backproject!(im, radonim) + im +end + +function backproject!(im::AbstractArray{T, 2}, radonim::AbstractArray{T, 2}) where T<:AbstractFloat + Nφ, Ns = size(radonim) + Nx, Ny = size(im) + Ldiag = hypot(Nx, Ny) + K = round(Int, Ldiag) + + dLz = Ldiag/(K-1) + dLs = Ldiag/(Ns-1) + dφ = π/Nφ + + im .= 0. + @inbounds for is = 1:Ns, iφ = 1:Nφ + φ = (iφ-1)*dφ + s = -Ldiag/2 + (is-1)*dLs + tx = cos(φ) + ty = sin(φ) + + for k=1:K + z = -Ldiag/2 + (k-1)*dLz + x = round(Int, Nx/2 + z*ty + s*tx) + y = round(Int, Ny/2 - z*tx + s*ty) + + if 1 <= x <= Nx && 1 <= y <= Ny + im[x, y] += radonim[iφ, is]*dLz + end + end + end + im +end + +end # module
