Mercurial > repos > PredictPDPS / changeset
changeset
merge
Thu, 25 Apr 2024 14:20:38 -0500
- author
- Tuomo Valkonen <tuomov@iki.fi>
- date
- Thu, 25 Apr 2024 14:20:38 -0500
- changeset 39
- b180a4f3b9bd
- parent 37
- bba159cf1636 (diff)
- parent 38
- 75116ad1d2e6 (current diff)
- child 40
- 2464329e356d
merge
| src/ImGenerate.jl | file | annotate | diff | comparison | revisions | |
| src/PET/PET.jl | file | annotate | diff | comparison | revisions |
--- a/Manifest.toml Thu Apr 25 19:16:17 2024 +0300 +++ b/Manifest.toml Thu Apr 25 14:20:38 2024 -0500 @@ -25,7 +25,7 @@ deps = ["DelimitedFiles", "Printf"] path = "../AlgTools" uuid = "c46e2e78-5339-41fd-a966-983ff60ab8e7" -version = "0.1.0" +version = "0.1.1" [[ArgTools]] uuid = "0dad84c5-d112-42e6-8d28-ef12dabb789f" @@ -40,12 +40,6 @@ [[Artifacts]] uuid = "56f22d72-fd6d-98f1-02f0-08ddc0907c33" -[[AssetRegistry]] -deps = ["Distributed", "JSON", "Pidfile", "SHA", "Test"] -git-tree-sha1 = "b25e88db7944f98789130d7b503276bc34bc098e" -uuid = "bf4720bc-e11a-5d0c-854e-bdca1663c893" -version = "0.1.0" - [[AxisAlgorithms]] deps = ["LinearAlgebra", "Random", "SparseArrays", "WoodburyMatrices"] git-tree-sha1 = "66771c8d21c8ff5e3a93379480a2307ac36863f7" @@ -66,12 +60,6 @@ uuid = "d1d4a3ce-64b1-5f1a-9ba4-7e7e69966f35" version = "0.1.8" -[[Blink]] -deps = ["Base64", "Distributed", "HTTP", "JSExpr", "JSON", "Lazy", "Logging", "MacroTools", "Mustache", "Mux", "Pkg", "Reexport", "Sockets", "WebIO"] -git-tree-sha1 = "bc93511973d1f949d45b0ea17878e6cb0ad484a1" -uuid = "ad839575-38b3-5650-b840-f874b8c74a25" -version = "0.12.9" - [[BufferedStreams]] git-tree-sha1 = "4ae47f9a4b1dc19897d3743ff13685925c5202ec" uuid = "e1450e63-4bb3-523b-b2a4-4ffa8c0fd77d" @@ -287,12 +275,6 @@ uuid = "2e619515-83b5-522b-bb60-26c02a35a201" version = "2.5.0+0" -[[EzXML]] -deps = ["Printf", "XML2_jll"] -git-tree-sha1 = "380053d61bb9064d6aa4a9777413b40429c79901" -uuid = "8f5d6c58-4d21-5cfd-889c-e3ad7ee6a615" -version = "1.2.0" - [[FFMPEG_jll]] deps = ["Artifacts", "Bzip2_jll", "FreeType2_jll", "FriBidi_jll", "JLLWrappers", "LAME_jll", "Libdl", "Ogg_jll", "OpenSSL_jll", "Opus_jll", "PCRE2_jll", "Zlib_jll", "libaom_jll", "libass_jll", "libfdk_aac_jll", "libvorbis_jll", "x264_jll", "x265_jll"] git-tree-sha1 = "ab3f7e1819dba9434a3a5126510c8fda3a4e7000" @@ -356,12 +338,6 @@ uuid = "559328eb-81f9-559d-9380-de523a88c83c" version = "1.0.10+0" -[[FunctionalCollections]] -deps = ["Test"] -git-tree-sha1 = "04cb9cfaa6ba5311973994fe3496ddec19b6292a" -uuid = "de31a74c-ac4f-5751-b3fd-e18cd04993ca" -version = "0.5.0" - [[Future]] deps = ["Random"] uuid = "9fa8497b-333b-5362-9e8d-4d0656e87820" @@ -450,12 +426,6 @@ uuid = "2e76f6c2-a576-52d4-95c1-20adfe4de566" version = "2.8.1+1" -[[Hiccup]] -deps = ["MacroTools", "Test"] -git-tree-sha1 = "6187bb2d5fcbb2007c39e7ac53308b0d371124bd" -uuid = "9fb69e20-1954-56bb-a84f-559cc56a8ff7" -version = "0.2.2" - [[Hwloc_jll]] deps = ["Artifacts", "JLLWrappers", "Libdl"] git-tree-sha1 = "ca0f6bf568b4bfc807e7537f081c81e35ceca114" @@ -547,10 +517,10 @@ version = "0.3.8" [[ImageTools]] -deps = ["AlgTools", "ColorTypes", "FileIO", "GR", "OffsetArrays", "Printf", "Setfield"] +deps = ["AlgTools", "ColorTypes", "ColorVectorSpace", "FileIO", "GR", "OffsetArrays", "Printf", "Setfield", "TestImages"] path = "../ImageTools" uuid = "b548cc0d-4ade-417e-bf62-0e39f9d2eee9" -version = "0.1.0" +version = "0.1.1" [[ImageTransformations]] deps = ["AxisAlgorithms", "ColorVectorSpace", "CoordinateTransformations", "ImageBase", "ImageCore", "Interpolations", "OffsetArrays", "Rotations", "StaticArrays"] @@ -655,12 +625,6 @@ uuid = "692b3bcd-3c85-4b1f-b108-f13ce0eb3210" version = "1.5.0" -[[JSExpr]] -deps = ["JSON", "MacroTools", "Observables", "WebIO"] -git-tree-sha1 = "b413a73785b98474d8af24fd4c8a975e31df3658" -uuid = "97c1335a-c9c5-57fe-bc5d-ec35cebe8660" -version = "0.5.4" - [[JSON]] deps = ["Dates", "Mmap", "Parsers", "Unicode"] git-tree-sha1 = "31e996f0a15c7b280ba9f76636b3ff9e2ae58c9a" @@ -679,12 +643,6 @@ uuid = "aacddb02-875f-59d6-b918-886e6ef4fbf8" version = "3.0.2+0" -[[Kaleido_jll]] -deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"] -git-tree-sha1 = "43032da5832754f58d14a91ffbe86d5f176acda9" -uuid = "f7e6163d-2fa5-5f23-b69c-1db539e41963" -version = "0.2.1+0" - [[LAME_jll]] deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"] git-tree-sha1 = "f6250b16881adf048549549fba48b1161acdac8c" @@ -714,12 +672,6 @@ uuid = "b964fa9f-0449-5b57-a5c2-d3ea65f4040f" version = "1.3.1" -[[Lazy]] -deps = ["MacroTools"] -git-tree-sha1 = "1370f8202dac30758f3c345f9909b97f53d87d3f" -uuid = "50d2b5c4-7a5e-59d5-8109-a42b560f39c0" -version = "0.15.1" - [[LazyArtifacts]] deps = ["Artifacts", "Pkg"] uuid = "4af54fe1-eca0-43a8-85a7-787d91b784e3" @@ -732,12 +684,12 @@ [[LibCURL]] deps = ["LibCURL_jll", "MozillaCACerts_jll"] uuid = "b27032c2-a3e7-50c8-80cd-2d36dbcbfd21" -version = "0.6.4" +version = "0.6.3" [[LibCURL_jll]] deps = ["Artifacts", "LibSSH2_jll", "Libdl", "MbedTLS_jll", "Zlib_jll", "nghttp2_jll"] uuid = "deac9b47-8bc7-5906-a0fe-35ac56dc84c0" -version = "8.4.0+0" +version = "7.84.0+0" [[LibGit2]] deps = ["Base64", "NetworkOptions", "Printf", "SHA"] @@ -746,7 +698,7 @@ [[LibSSH2_jll]] deps = ["Artifacts", "Libdl", "MbedTLS_jll"] uuid = "29816b5a-b9ab-546f-933c-edad1886dfa8" -version = "1.11.0+1" +version = "1.10.2+0" [[Libdl]] uuid = "8f399da3-3557-5675-b5ff-fb832c97cbdb" @@ -921,18 +873,6 @@ uuid = "14a3606d-f60d-562e-9121-12d972cd8159" version = "2022.10.11" -[[Mustache]] -deps = ["Printf", "Tables"] -git-tree-sha1 = "a7cefa21a2ff993bff0456bf7521f46fc077ddf1" -uuid = "ffc61752-8dc7-55ee-8c37-f3e9cdd09e70" -version = "1.0.19" - -[[Mux]] -deps = ["AssetRegistry", "Base64", "HTTP", "Hiccup", "MbedTLS", "Pkg", "Sockets"] -git-tree-sha1 = "7295d849103ac4fcbe3b2e439f229c5cc77b9b69" -uuid = "a975b10e-0019-58db-a62f-e48ff68538c9" -version = "1.0.2" - [[NaNMath]] deps = ["OpenLibm_jll"] git-tree-sha1 = "0877504529a3e5c3343c6f8b4c0381e57e4387e4" @@ -955,11 +895,6 @@ uuid = "ca575930-c2e3-43a9-ace4-1e988b2c1908" version = "1.2.0" -[[Observables]] -git-tree-sha1 = "7438a59546cf62428fc9d1bc94729146d37a7225" -uuid = "510215fc-4207-5dde-b226-833fc4488ee2" -version = "0.5.5" - [[OffsetArrays]] git-tree-sha1 = "e64b4f5ea6b7389f6f046d13d4896a8f9c1ba71e" uuid = "6fe1bfb0-de20-5000-8ca7-80f57d26f881" @@ -1073,12 +1008,6 @@ uuid = "54e51dfa-9dd7-5231-aa84-a4037b83483a" version = "0.3.6" -[[Pidfile]] -deps = ["FileWatching", "Test"] -git-tree-sha1 = "2d8aaf8ee10df53d0dfb9b8ee44ae7c04ced2b03" -uuid = "fa939f87-e72e-5be4-a000-7fc836dbe307" -version = "1.3.0" - [[Pixman_jll]] deps = ["Artifacts", "CompilerSupportLibraries_jll", "JLLWrappers", "LLVMOpenMP_jll", "Libdl"] git-tree-sha1 = "64779bc4c9784fee475689a1752ef4d5747c5e87" @@ -1096,36 +1025,6 @@ uuid = "eebad327-c553-4316-9ea0-9fa01ccd7688" version = "0.3.3" -[[PlotlyBase]] -deps = ["ColorSchemes", "Dates", "DelimitedFiles", "DocStringExtensions", "JSON", "LaTeXStrings", "Logging", "Parameters", "Pkg", "REPL", "Requires", "Statistics", "UUIDs"] -git-tree-sha1 = "56baf69781fc5e61607c3e46227ab17f7040ffa2" -uuid = "a03496cd-edff-5a9b-9e67-9cda94a718b5" -version = "0.8.19" - -[[PlotlyJS]] -deps = ["Base64", "Blink", "DelimitedFiles", "JSExpr", "JSON", "Kaleido_jll", "Markdown", "Pkg", "PlotlyBase", "PlotlyKaleido", "REPL", "Reexport", "Requires", "WebIO"] -git-tree-sha1 = "e62d886d33b81c371c9d4e2f70663c0637f19459" -uuid = "f0f68f2c-4968-5e81-91da-67840de0976a" -version = "0.18.13" - - [PlotlyJS.extensions] - CSVExt = "CSV" - DataFramesExt = ["DataFrames", "CSV"] - IJuliaExt = "IJulia" - JSON3Ext = "JSON3" - - [PlotlyJS.weakdeps] - CSV = "336ed68f-0bac-5ca0-87d4-7b16caf5d00b" - DataFrames = "a93c6f00-e57d-5684-b7b6-d8193f3e46c0" - IJulia = "7073ff75-c697-5162-941a-fcdaad2a7d2a" - JSON3 = "0f8b85d8-7281-11e9-16c2-39a750bddbf1" - -[[PlotlyKaleido]] -deps = ["Base64", "JSON", "Kaleido_jll"] -git-tree-sha1 = "2650cd8fb83f73394996d507b3411a7316f6f184" -uuid = "f2990250-8cf9-495f-b13a-cce12b45703c" -version = "2.2.4" - [[PoissonRandom]] deps = ["Random"] git-tree-sha1 = "a0f1159c33f846aa77c3f30ebbc69795e5327152" @@ -1496,24 +1395,6 @@ uuid = "ea10d353-3f73-51f8-a26c-33c1cb351aa5" version = "1.4.2" -[[WebIO]] -deps = ["AssetRegistry", "Base64", "Distributed", "FunctionalCollections", "JSON", "Logging", "Observables", "Pkg", "Random", "Requires", "Sockets", "UUIDs", "WebSockets", "Widgets"] -git-tree-sha1 = "0eef0765186f7452e52236fa42ca8c9b3c11c6e3" -uuid = "0f1e0344-ec1d-5b48-a673-e5cf874b6c29" -version = "0.8.21" - -[[WebSockets]] -deps = ["Base64", "Dates", "HTTP", "Logging", "Sockets"] -git-tree-sha1 = "4162e95e05e79922e44b9952ccbc262832e4ad07" -uuid = "104b5d7c-a370-577a-8038-80a2059c5097" -version = "1.6.0" - -[[Widgets]] -deps = ["Colors", "Dates", "Observables", "OrderedCollections"] -git-tree-sha1 = "fcdae142c1cfc7d89de2d11e08721d0f2f86c98a" -uuid = "cc8bc4a8-27d6-5769-a93b-9d913e69aa62" -version = "0.6.6" - [[WoodburyMatrices]] deps = ["LinearAlgebra", "SparseArrays"] git-tree-sha1 = "5f24e158cf4cee437052371455fe361f526da062" @@ -1525,12 +1406,6 @@ uuid = "76eceee3-57b5-4d4a-8e66-0e911cebbf60" version = "1.6.1" -[[XLSX]] -deps = ["Artifacts", "Dates", "EzXML", "Printf", "Tables", "ZipFile"] -git-tree-sha1 = "319b05e790046f18f12b8eae542546518ef1a88f" -uuid = "fdbf4ff8-1666-58a4-91e7-1b58723a45e0" -version = "0.10.1" - [[XML2_jll]] deps = ["Artifacts", "JLLWrappers", "Libdl", "Libiconv_jll", "Zlib_jll"] git-tree-sha1 = "532e22cf7be8462035d092ff21fada7527e2c488" @@ -1693,12 +1568,6 @@ uuid = "c5fb5394-a638-5e4d-96e5-b29de1b5cf10" version = "1.5.0+0" -[[ZipFile]] -deps = ["Libdl", "Printf", "Zlib_jll"] -git-tree-sha1 = "f492b7fe1698e623024e873244f10d89c95c340a" -uuid = "a5390f91-8eb1-5f08-bee0-b1d1ffed6cea" -version = "0.10.1" - [[Zlib_jll]] deps = ["Libdl"] uuid = "83775a58-1f1d-513f-b197-d71354ab007a" @@ -1790,7 +1659,7 @@ [[nghttp2_jll]] deps = ["Artifacts", "Libdl"] uuid = "8e850ede-7688-5339-a07c-302acd2aaf8d" -version = "1.52.0+1" +version = "1.48.0+0" [[oneTBB_jll]] deps = ["Artifacts", "JLLWrappers", "Libdl"]
--- a/Project.toml Thu Apr 25 19:16:17 2024 +0300 +++ b/Project.toml Thu Apr 25 14:20:38 2024 -0500 @@ -22,7 +22,6 @@ MAT = "23992714-dd62-5051-b70f-ba57cb901cac" OffsetArrays = "6fe1bfb0-de20-5000-8ca7-80f57d26f881" PerceptualColourMaps = "54e51dfa-9dd7-5231-aa84-a4037b83483a" -PlotlyJS = "f0f68f2c-4968-5e81-91da-67840de0976a" PoissonRandom = "e409e4f3-bfea-5376-8464-e040bb5c01ab" Printf = "de0858da-6303-5e67-8744-51eddeeeb8d7" QuartzImageIO = "dca85d43-d64c-5e67-8c65-017450d5d020" @@ -32,4 +31,3 @@ StableRNGs = "860ef19b-820b-49d6-a774-d7a799459cd3" Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2" TestImages = "5e47fb64-e119-507b-a336-dd2b206d9990" -XLSX = "fdbf4ff8-1666-58a4-91e7-1b58723a45e0"
--- a/src/Algorithm.jl Thu Apr 25 19:16:17 2024 +0300 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,148 +0,0 @@ -#################################################################### -# Predictive online PDPS for optical flow with known velocity field -#################################################################### - -__precompile__() - -module Algorithm - -identifier = "pdps_known" - -using Printf - -using AlgTools.Util -import AlgTools.Iterate -using ImageTools.Gradient - -using ..OpticalFlow: ImageSize, - Image, - pdflow! - -######################### -# Iterate initialisation -######################### - -function init_rest(x::Image) - imdim=size(x) - - y = zeros(2, imdim...) - Δx = copy(x) - Δy = copy(y) - x̄ = copy(x) - - return x, y, Δx, Δy, x̄ -end - -function init_iterates(xinit::Image) - return init_rest(copy(xinit)) -end - -function init_iterates(dim::ImageSize) - return init_rest(zeros(dim...)) -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+γ*τ ≥ Λ) - σ = σ₀*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²) - - ###################### - # Initialise iterates - ###################### - - x, y, Δx, Δy, x̄ = init_iterates(dim) - init_data = (params.init == :data) - - #################### - # Run the algorithm - #################### - - v = iterate(params) do verbose :: Function, - b :: Image, - v_known :: DisplacementT, - 🚫unused_b_next :: Image - - ################## - # Prediction step - ################## - if init_data - x .= b - init_data = false - end - - pdflow!(x, Δx, y, Δy, v_known, params.dual_flow) - - if params.prox_predict - ∇₂!(Δy, x) - @. y = (y + σ̃*Δy)/(1 + σ̃*(ρ̃+ρ/α)) - proj_norm₂₁ball!(y, α) - end - - ############ - # PDPS step - ############ - - ∇₂ᵀ!(Δx, y) # primal step: - @. x̄ = x # | save old x for over-relax - @. x = (x-τ*(Δx-b))/(1+τ) # | prox - @. x̄ = 2x - x̄ # over-relax - ∇₂!(Δy, x̄) # dual step: y - @. y = (y + σ*Δy)/(1 + σ*ρ/α) # | - proj_norm₂₁ball!(y, α) # | prox - - ################################ - # Give function value if needed - ################################ - v = verbose() do - ∇₂!(Δy, x) - value = norm₂²(b-x)/2 + params.α*γnorm₂₁(Δy, params.ρ) - value, x, [NaN, NaN], nothing - end - - v - end - - return x, y, v -end - -end # Module - -
--- a/src/AlgorithmBoth.jl Thu Apr 25 19:16:17 2024 +0300 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,239 +0,0 @@ -###################################################################### -# Predictive online PDPS for optical flow with unknown velocity field -###################################################################### - -__precompile__() - -module AlgorithmBoth - -identifier = "pdps_unknown_basic" - -using Printf - -using AlgTools.Util -import AlgTools.Iterate -using ImageTools.Gradient - -using ..OpticalFlow: ImageSize, - Image, - Gradient, - DisplacementConstant, - DisplacementFull, - pdflow!, - pointwise_gradiprod_2d!, - pointwise_gradiprod_2dᵀ!, - filter_hs - -using ..Algorithm: step_lengths - -############# -# Data types -############# - -struct Primal{DisplacementT} - x :: Image - u :: DisplacementT -end - -function Base.similar(x::Primal{DisplacementT}) where DisplacementT - return Primal{DisplacementT}(Base.similar(x.x), Base.similar(x.u)) -end - -function Base.copy(x::Primal{DisplacementT}) where DisplacementT - return Primal{DisplacementT}(Base.copy(x.x), Base.copy(x.u)) - end - -struct Dual - tv :: Gradient - flow :: Image -end - -function Base.similar(y::Dual) - return Dual(Base.similar(y.tv), Base.similar(y.flow)) -end - -function Base.copy(y::Dual) - return Dual(Base.copy(y.tv), Base.copy(y.flow)) - end - -######################### -# Iterate initialisation -######################### - -function init_primal(xinit::Image, ::Type{DisplacementConstant}) - return Primal{DisplacementConstant}(xinit, zeros(2)) -end - -function init_primal(xinit::Image, ::Type{DisplacementFull}) - return Primal{DisplacementFull}(xinit, zeros(2, size(xinit)...)) -end - -function init_rest(x::Primal{DisplacementT}) where DisplacementT - imdim=size(x.x) - - y = Dual(zeros(2, imdim...), zeros(imdim)) - Δx = copy(x) - Δy = copy(y) - x̄ = copy(x) - - return x, y, Δx, Δy, x̄ -end - -function init_iterates( :: Type{DisplacementT}, - xinit::Primal{DisplacementT}) where DisplacementT - return init_rest(copy(xinit)) -end - -function init_iterates( :: Type{DisplacementT}, xinit::Image) where DisplacementT - return init_rest(init_primal(copy(xinit), DisplacementT)) -end - -function init_iterates( :: Type{DisplacementT}, dim::ImageSize) where DisplacementT - return init_rest(init_primal(zeros(dim...), DisplacementT)) -end - -############################################## -# Weighting for different displacements types -############################################## - -norm²weight( :: Type{DisplacementConstant}, sz ) = prod(sz) -norm²weight( :: Type{DisplacementFull}, sz ) = 1 - -############ -# Algorithm -############ - -function solve( :: Type{DisplacementT}; - dim :: ImageSize, - iterate = AlgTools.simple_iterate, - params::NamedTuple) where DisplacementT - - ###################### - # Initialise iterates - ###################### - - x, y, Δx, Δy, x̄ = init_iterates(DisplacementT, dim) - init_data = (params.init == :data) - - # … for tracking cumulative movement - if DisplacementT == DisplacementConstant - ucumul = [0.0, 0.0] - else - ucumul = [NaN, NaN] - end - - ############################################# - # Extract parameters and set up step lengths - ############################################# - - α, ρ, λ, θ, T = params.α, params.ρ, params.λ, params.θ, params.timestep - R_K² = max(∇₂_norm₂₂_est², ∇₂_norm₂∞_est²*params.dynrange^2) - γ = min(1, λ*norm²weight(DisplacementT, size(x.x))) - τ, σ, σ̃, ρ̃ = step_lengths(params, γ, R_K²) - - kernel = params.kernel - - #################### - # Run the algorithm - #################### - - b_next_filt=nothing - - v = iterate(params) do verbose :: Function, - b :: Image, - 🚫unused_v_known :: DisplacementT, - b_next :: Image - - #################################### - # Smooth data for Horn–Schunck term - #################################### - - b_filt, b_next_filt = filter_hs(b, b_next, b_next_filt, kernel) - - ############################ - # Construct K for this step - ############################ - - K! = (yʹ, xʹ) -> begin - # Optical flow part: ⟨⟨u, ∇b_k⟩⟩. - # Use y.tv as temporary gradient storage. - pointwise_gradiprod_2d!(yʹ.flow, yʹ.tv, xʹ.u, b_filt) - #@. yʹ.flow = -yʹ.flow - # TV part - ∇₂!(yʹ.tv, xʹ.x) - end - Kᵀ! = (xʹ, yʹ) -> begin - # Optical flow part: ∇b_k ⋅ y - pointwise_gradiprod_2dᵀ!(xʹ.u, yʹ.flow, b_filt) - #@. xʹ.u = -xʹ.u - # TV part - ∇₂ᵀ!(xʹ.x, yʹ.tv) - end - - ################## - # Prediction step - ################## - - if init_data - x .= b - init_data = false - end - - pdflow!(x.x, Δx.x, y.tv, Δy.tv, y.flow, Δy.flow, x.u, params.dual_flow) - - # Predict zero displacement - x.u .= 0 - if params.prox_predict - K!(Δy, x) - @. y.tv = (y.tv + σ̃*Δy.tv)/(1 + σ̃*(ρ̃+ρ/α)) - proj_norm₂₁ball!(y.tv, α) - @. y.flow = (y.flow+σ̃*((b_next_filt-b_filt)/T+Δy.flow))/(1+σ̃*(ρ̃+1/θ)) - end - - ############ - # PDPS step - # - # NOTE: For DisplacementConstant, the x.u update is supposed to be with - # respect to the 𝟙^*𝟙 norm/inner product that makes the norm equivalent - # to full-space norm when restricted to constant displacements. Since - # `OpticalFlow.pointwise_gradiprod_2dᵀ!` already uses this inner product, - # and the λ-weighted term in the problem is with respect to this norm, - # all the norm weights disappear in this update. - ############ - - Kᵀ!(Δx, y) # primal step: - @. x̄.x = x.x # | save old x for over-relax - @. x̄.u = x.u # | - @. x.x = (x.x-τ*(Δx.x-b))/(1+τ) # | prox - @. x.u = (x.u-τ*Δx.u)/(1+τ*λ) # | - @. x̄.x = 2x.x - x̄.x # over-relax - @. x̄.u = 2x.u - x̄.u # | - K!(Δy, x̄) # dual step: y - @. y.tv = (y.tv + σ*Δy.tv)/(1 + σ*ρ/α) # | - proj_norm₂₁ball!(y.tv, α) # | prox - @. y.flow = (y.flow+σ*((b_next_filt-b_filt)/T+Δy.flow))/(1+σ/θ) - - if DisplacementT == DisplacementConstant - ucumul .+= x.u - end - - ######################################################## - # Give function value and cumulative movement if needed - ######################################################## - v = verbose() do - K!(Δy, x) - value = (norm₂²(b-x.x)/2 + θ*norm₂²((b_next_filt-b_filt)./T+Δy.flow) - + λ*norm₂²(x.u)/2 + α*γnorm₂₁(Δy.tv, ρ)) - - value, x.x, ucumul, nothing - end - - return v - end - - return x, y, v -end - -end # Module - -
--- a/src/AlgorithmBothCumul.jl Thu Apr 25 19:16:17 2024 +0300 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,139 +0,0 @@ -###################################################################### -# Predictive online PDPS for optical flow with unknown velocity field -###################################################################### - -__precompile__() - -module AlgorithmBothCumul - -identifier = "pdps_unknown_cumul" - -using Printf - -using AlgTools.Util -import AlgTools.Iterate -using ImageTools.Gradient -using ImageTools.ImFilter - -using ..OpticalFlow: Image, - ImageSize, - DisplacementConstant, - pdflow!, - horn_schunck_reg_prox!, - pointwise_gradiprod_2d! - -using ..AlgorithmBothGreedyV: init_iterates -using ..Algorithm: step_lengths - -############ -# Algorithm -############ - -function solve( :: Type{DisplacementT}; - dim :: ImageSize, - iterate = AlgTools.simple_iterate, - params::NamedTuple) where DisplacementT - - ###################### - # Initialise iterates - ###################### - - x, y, Δx, Δy, x̄, u = init_iterates(DisplacementT, dim) - init_data = (params.init == :data) - - # … for tracking cumulative movement - if DisplacementT == DisplacementConstant - ucumul = zeros(size(u)...) - end - - ############################################# - # Extract parameters and set up step lengths - ############################################# - - α, ρ, λ, θ, T = params.α, params.ρ, params.λ, params.θ, params.timestep - R_K² = ∇₂_norm₂₂_est² - γ = 1 - τ, σ, σ̃, ρ̃ = step_lengths(params, γ, R_K²) - - kernel = params.kernel - - #################### - # Run the algorithm - #################### - - b₀=nothing - b₀_filt=nothing - u_prev=zeros(size(u)) - - v = iterate(params) do verbose :: Function, - b :: Image, - 🚫unused_v_known :: DisplacementT, - 🚫unused_b_next :: Image - - ######################################################### - # Smoothen data for Horn–Schunck term; zero initial data - ######################################################### - - b_filt = (kernel==nothing ? b : simple_imfilter(b, kernel)) - - if b₀ == nothing - b₀ = b - b₀_filt = b_filt - end - - ################################################ - # Prediction step - # We leave u as-is in this cumulative version - ################################################ - - if init_data - x .= b - init_data = false - end - - pdflow!(x, Δx, y, Δy, u-u_prev, params.dual_flow) - - if params.prox_predict - ∇₂!(Δy, x) - @. y = (y + σ̃*Δy)/(1 + σ̃*(ρ̃+ρ/α)) - proj_norm₂₁ball!(y, α) - end - - # Store current cumulative displacement before updating in next step. - u_prev .= u - - ############ - # PDPS step - ############ - - ∇₂ᵀ!(Δx, y) # primal step: - @. x̄ = x # | save old x for over-relax - @. x = (x-τ*(Δx-b))/(1+τ) # | prox - horn_schunck_reg_prox!(u, b_filt, b₀_filt, τ, θ, λ, T) - @. x̄ = 2x - x̄ # over-relax - ∇₂!(Δy, x̄) # dual step: y - @. y = (y + σ*Δy)/(1 + σ*ρ/α) # | - proj_norm₂₁ball!(y, α) # | prox - - ######################################################## - # Give function value and cumulative movement if needed - ######################################################## - v = verbose() do - ∇₂!(Δy, x) - tmp = zeros(size(b_filt)) - pointwise_gradiprod_2d!(tmp, Δy, u, b₀_filt) - value = (norm₂²(b-x)/2 + θ*norm₂²((b_filt-b₀_filt)./T+tmp) - + λ*norm₂²(u)/2 + α*γnorm₂₁(Δy, ρ)) - - value, x, u, nothing - end - - return v - end - - return x, y, v -end - -end # Module - -
--- a/src/AlgorithmBothGreedyV.jl Thu Apr 25 19:16:17 2024 +0300 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,167 +0,0 @@ -###################################################################### -# Predictive online PDPS for optical flow with unknown velocity field -###################################################################### - -__precompile__() - -module AlgorithmBothGreedyV - -identifier = "pdps_unknown_greedyv" - -using Printf - -using AlgTools.Util -import AlgTools.Iterate -using ImageTools.Gradient - -using ..OpticalFlow: Image, - ImageSize, - DisplacementConstant, - DisplacementFull, - pdflow!, - horn_schunck_reg_prox!, - pointwise_gradiprod_2d!, - filter_hs - -using ..Algorithm: step_lengths - -######################### -# Iterate initialisation -######################### - -function init_displ(xinit::Image, ::Type{DisplacementConstant}) - return xinit, zeros(2) -end - -function init_displ(xinit::Image, ::Type{DisplacementFull}) - return xinit, zeros(2, size(xinit)...) -end - -function init_rest(x::Image, u::DisplacementT) where DisplacementT - imdim=size(x) - - y = zeros(2, imdim...) - Δx = copy(x) - Δy = copy(y) - x̄ = copy(x) - - return x, y, Δx, Δy, x̄, u -end - -function init_iterates( :: Type{DisplacementT}, xinit::Image) where DisplacementT - return init_rest(init_displ(copy(xinit), DisplacementT)...) -end - -function init_iterates( :: Type{DisplacementT}, dim::ImageSize) where DisplacementT - return init_rest(init_displ(zeros(dim...), DisplacementT)...) -end - -############ -# Algorithm -############ - -function solve( :: Type{DisplacementT}; - dim :: ImageSize, - iterate = AlgTools.simple_iterate, - params::NamedTuple) where DisplacementT - - ###################### - # Initialise iterates - ###################### - - x, y, Δx, Δy, x̄, u = init_iterates(DisplacementT, dim) - init_data = (params.init == :data) - - # … for tracking cumulative movement - if DisplacementT == DisplacementConstant - ucumul = [0.0, 0.0] - else - ucumul = [NaN, NaN] - end - - ############################################# - # Extract parameters and set up step lengths - ############################################# - - α, ρ, λ, θ, T = params.α, params.ρ, params.λ, params.θ, params.timestep - R_K² = ∇₂_norm₂₂_est² - γ = 1 - τ, σ, σ̃, ρ̃ = step_lengths(params, γ, R_K²) - - kernel = params.kernel - - #################### - # Run the algorithm - #################### - - b_next_filt=nothing - - v = iterate(params) do verbose :: Function, - b :: Image, - 🚫unused_v_known :: DisplacementT, - b_next :: Image - - #################################### - # Smooth data for Horn–Schunck term - #################################### - - b_filt, b_next_filt = filter_hs(b, b_next, b_next_filt, kernel) - - ################## - # Prediction step - ################## - - if init_data - x .= b - init_data = false - end - - pdflow!(x, Δx, y, Δy, u, params.dual_flow) - - # Predict zero displacement - u .= 0 - if params.prox_predict - ∇₂!(y, x) - @. y = (y + σ̃*Δy)/(1 + σ̃*(ρ̃+ρ/α)) - proj_norm₂₁ball!(y, α) - end - - ############ - # PDPS step - ############ - - ∇₂ᵀ!(Δx, y) # primal step: - @. x̄ = x # | save old x for over-relax - @. x = (x-τ*(Δx-b))/(1+τ) # | prox - horn_schunck_reg_prox!(u, b_next_filt, b_filt, θ, λ, T, τ) - @. x̄ = 2x - x̄ # over-relax - ∇₂!(y, x̄) # dual step: y - @. y = (y + σ*Δy)/(1 + σ*ρ/α) # | - proj_norm₂₁ball!(y, α) # | prox - - if DisplacementT == DisplacementConstant - ucumul .+= u - end - - ######################################################## - # Give function value and cumulative movement if needed - ######################################################## - v = verbose() do - ∇₂!(Δy, x) - tmp = zeros(size(b_filt)) - pointwise_gradiprod_2d!(tmp, Δy, u, b_filt) - value = (norm₂²(b-x)/2 + θ*norm₂²((b_next_filt-b_filt)./T+tmp) - + λ*norm₂²(u)/2 + α*γnorm₂₁(Δy, ρ)) - - value, x, ucumul, nothing - end - - return v - end - - return x, y, v -end - -end # Module - -
--- a/src/AlgorithmBothMulti.jl Thu Apr 25 19:16:17 2024 +0300 +++ b/src/AlgorithmBothMulti.jl Thu Apr 25 14:20:38 2024 -0500 @@ -26,7 +26,7 @@ ConstantDisplacementHornSchunckData, filter_hs -using ..Algorithm: step_lengths +using ..AlgorithmProximal: step_lengths ######################### # Iterate initialisation
--- a/src/AlgorithmBothMulti.jl.orig Thu Apr 25 19:16:17 2024 +0300 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,277 +0,0 @@ -###################################################################### -# Predictive online PDPS for optical flow with unknown velocity field -###################################################################### - -__precompile__() - -module AlgorithmBothMulti - -using Printf - -using AlgTools.Util -import AlgTools.Iterate -using ImageTools.Gradient -using ImageTools.ImFilter - -using ..OpticalFlow: Image, - ImageSize, - DisplacementConstant, - pdflow!, - horn_schunck_reg_prox!, - pointwise_gradiprod_2d!, - horn_schunck_reg_prox_op!, - mldivide_step_plus_sym2x2!, - ConstantDisplacementHornSchunckData, - filter_hs - -using ..Algorithm: step_lengths - -######################### -# Iterate initialisation -######################### - -function init_displ(xinit::Image, ::Type{DisplacementConstant}, n::Integer) - return xinit, zeros(n, 2) -end - -# function init_displ(xinit::Image, ::Type{DisplacementFull}, n::Integer) -# return xinit, zeros(n, 2, size(xinit)...) -# end - -function init_rest(x::Image, u::DisplacementT) where DisplacementT - imdim=size(x) - - y = zeros(2, imdim...) - Δx = copy(x) - Δy = copy(y) - x̄ = copy(x) - - return x, y, Δx, Δy, x̄, u -end - -function init_iterates( :: Type{DisplacementT}, - xinit::Image, - n::Integer) where DisplacementT - return init_rest(init_displ(copy(xinit), DisplacementT, n)...) -end - -function init_iterates( :: Type{DisplacementT}, - dim::ImageSize, - n::Integer) where DisplacementT - return init_rest(init_displ(zeros(dim...), DisplacementT, n)...) -end - -############ -# Algorithm -############ - -function solve( :: Type{DisplacementT}; - dim :: ImageSize, - iterate = AlgTools.simple_iterate, - params::NamedTuple) where DisplacementT - - ###################### - # Initialise iterates - ###################### - - n = params.displacement_count - k = 0 # number of displacements we have already - - x, y, Δx, Δy, x̄, u = init_iterates(DisplacementT, dim, n) - init_data = (params.init == :data) - hs = [ConstantDisplacementHornSchunckData() for i=1:n] - #hs = Array{ConstantDisplacementHornSchunckData}(undef, n) - A = Array{Float64,3}(undef, n, 2, 2) - d = Array{Float64,2}(undef, n, 2) - - # … for tracking cumulative movement - ucumulbase = [0.0, 0.0] - - ############################################# - # Extract parameters and set up step lengths - ############################################# - - α, ρ, λ, θ = params.α, params.ρ, params.λ, params.θ - R_K² = ∇₂_norm₂₂_est² - γ = 1 - τ, σ, σ̃, ρ̃ = step_lengths(params, γ, R_K²) - - kernel = params.kernel - T = params.timestep - - #################### - # Run the algorithm - #################### - - b_next_filt = nothing - diffu = similar(u[1, :]) - - v = iterate(params) do verbose :: Function, - b :: Image, - 🚫unused_v_known :: DisplacementT, - b_next :: Image - - #################################### - # Smooth data for Horn–Schunck term - #################################### - - b_filt, b_next_filt = filter_hs(b, b_next, b_next_filt, kernel) - - ################################################ - # Prediction step - ################################################ - - # Predict x and y - if k==0 - if init_data - x .= b - init_data = false - end - else - # Displacement from previous to this image is estimated as - # the difference of their displacements from beginning of window. - if k>1 - @. @views diffu = u[k, :] - u[k-1, :] - else - @. @views diffu = u[k, :] - end - - pdflow!(x, Δx, y, Δy, diffu, params.dual_flow) - end - - # Shift stored prox matrices - if k==n - tmp = copy(u[1, :]) - ucumulbase .+= tmp - for j=1:(n-1) - @. @views u[j, :] = u[j+1, :] - tmp - hs[j] = hs[j+1] - end - # Create new struct as original contains references to objects that - # have been moved to index n-1. - hs[n]=ConstantDisplacementHornSchunckData() - else - k += 1 - end - - # Predict u: zero displacement from current to next image, i.e., - # same displacement to beginning of window. - if k==1 - @. @views u[k, :] = 0.0 - else - @. @views u[k, :] = u[k-1, :] - end - - # Predictor proximals tep - if params.prox_predict - ∇₂!(Δy, x) - @. y = (y + σ̃*Δy)/(1 + σ̃*(ρ̃+ρ/α)) - proj_norm₂₁ball!(y, α) - end - - ################################################################################# - # PDPS step - # - # For the displacements, with τ̃=τ/k, we need to solve for 2≤j<k, - # - # (1) (I/τ̃+M₀^j+M₀^{j+1})u^j = M₀^ju^{j-1} + M₀^{j+1}u^{j+1} - # + ũ^j/τ̃ - z^j + z^{j+1}, - # - # as well as - # - # (2) (I/τ̃+M₀^k)u^k = M₀^k u^{k-1} + ũ^k/τ̃ - z^k - # - # and - # - # (3) (I/τ̃+M₀^1+M₀^2)u^1 = 0 + M₀^{2}u^{2} + ũ^1/τ̃ - z^1 + z^{2} - # - # We first construct from (2) that - # - # u^k = A^k u^{k-1} + d^k - # - # for - # - # A^k := (I/τ̃+M₀^k)^{-1} M₀^k - # d_k := (I/τ̃+M₀^k)^{-1} (ũ^k/τ̃ - z^k). - # - # Inserting this into (1) we need - # - # (4) (I/τ̃+M₀^j+M₀^{j+1}(I-A^{j+1}))u^j = M₀^ju^{j-1} + M₀^{j+1}d^{j+1} - # + ũ^j/τ̃ - z^j + z^{j+1}. - # - # This is well-defined because A^{j+1} < I. It also has the same form as (1), so - # we continue with - # - # (5) u^j = A^j u^{j-1} + d^j - # - # for - # - # A^j := (I/τ̃+M₀^j+M₀^{j+1}(I-A^{j+1}))^{-1} M₀^j - # d^j := (I/τ̃+M₀^j+M₀^{j+1}(I-A^{j+1}))^{-1} - # (M₀^{j+1}d^{j+1} + ũ^j/τ̃ - z^j + z^{j+1}) - # - # Finally from (3) with these we need - # - # (I/τ̃+M₀^1+M₀^2(I-A^2))u^1 = M₀^2d^2 + ũ^1/τ̃ - z^1 + z^2, - # - # which is of the same form as (4) with u^0=0, so by (5) u^1=d^1. - # - ################################################################################# - - ∇₂ᵀ!(Δx, y) # primal step: - @. x̄ = x # | save old x for over-relax - @. x = (x-τ*(Δx-b))/(1+τ) # | prox - # | | for displacement - # Calculate matrices for latest data; rest is stored. - @views begin - horn_schunck_reg_prox_op!(hs[k], b_next_filt, b_filt, θ, λ, T) - - τ̃=τ/k - - B = hs[k].M₀ - c = u[k, :]./τ̃-hs[k].z - mldivide_step_plus_sym2x2!(A[k, :, :], B, B, τ̃) - mldivide_step_plus_sym2x2!(d[k, :], B, c, τ̃) - - for j=(k-1):-1:1 - B = hs[j].M₀+hs[j+1].M₀*([1 0; 0 1]-A[j+1, :, :]) - c = hs[j+1].M₀*d[j+1, :]+u[j, :]./τ̃-hs[j].z+hs[j+1].z - mldivide_step_plus_sym2x2!(A[j, :, :], B, hs[j].M₀, τ̃) - mldivide_step_plus_sym2x2!(d[j, :], B, c, τ̃) - end - - u[1, :] .= d[1, :] - for j=2:k - u[j, :] .= A[j, :, :]*u[j-1, :] + d[j, :] - end - end - - @. x̄ = 2x - x̄ # over-relax: x̄ = 2x-x_old - ∇₂!(Δy, x̄) # dual step: y - @. y = (y + σ*Δy)/(1 + σ*ρ/α) # | - proj_norm₂₁ball!(y, α) # | prox - - ######################################################## - # Give function value and cumulative movement if needed - ######################################################## - v = verbose() do - ∇₂!(Δy, x) - hs_plus_reg=0 - for j=1:k - v=(j==1 ? u[j, :] : u[j, :]-u[j-1, :]) - hs_plus_reg += hs[j].cv/2 + dot(hs[j].Mv*v, v)/2+dot(hs[j].av, v) - end - value = (norm₂²(b-x)/2 + hs_plus_reg/k + α*γnorm₂₁(Δy, ρ)) - - value, x, u[k, :]+ucumulbase, u[1:k,:].+ucumulbase' - end - - return v - end - - return x, y, v -end - -end # Module - -
--- a/src/AlgorithmBothNL.jl Thu Apr 25 19:16:17 2024 +0300 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,242 +0,0 @@ -###################################################################### -# Predictive online PDPS for optical flow with unknown velocity field -###################################################################### - -__precompile__() - -module AlgorithmBothNL - -identifier = "pdps_unknown_nl" - - -using Printf - -using AlgTools.Util -import AlgTools.Iterate -using ImageTools.Gradient - -using ..OpticalFlow: ImageSize, - Image, - Gradient, - DisplacementConstant, - DisplacementFull, - pdflow!, - pointwise_gradiprod_2d!, - pointwise_gradiprod_2dᵀ!, - filter_hs - -using ..Algorithm: step_lengths - -############# -# Data types -############# - -struct Primal{DisplacementT} - x :: Image - u :: DisplacementT -end - -function Base.similar(x::Primal{DisplacementT}) where DisplacementT - return Primal{DisplacementT}(Base.similar(x.x), Base.similar(x.u)) -end - -function Base.copy(x::Primal{DisplacementT}) where DisplacementT - return Primal{DisplacementT}(Base.copy(x.x), Base.copy(x.u)) - end - -struct Dual - tv :: Gradient - flow :: Image -end - -function Base.similar(y::Dual) - return Dual(Base.similar(y.tv), Base.similar(y.flow)) -end - -function Base.copy(y::Dual) - return Dual(Base.copy(y.tv), Base.copy(y.flow)) - end - -######################### -# Iterate initialisation -######################### - -function init_primal(xinit::Image, ::Type{DisplacementConstant}) - return Primal{DisplacementConstant}(xinit, zeros(2)) -end - -function init_primal(xinit::Image, ::Type{DisplacementFull}) - return Primal{DisplacementFull}(xinit, zeros(2, size(xinit)...)) -end - -function init_rest(x::Primal{DisplacementT}) where DisplacementT - imdim=size(x.x) - - y = Dual(zeros(2, imdim...), zeros(imdim)) - Δx = copy(x) - Δy = copy(y) - x̄ = copy(x) - - return x, y, Δx, Δy, x̄ -end - -function init_iterates( :: Type{DisplacementT}, - xinit::Primal{DisplacementT}) where DisplacementT - return init_rest(copy(xinit)) -end - -function init_iterates( :: Type{DisplacementT}, xinit::Image) where DisplacementT - return init_rest(init_primal(copy(xinit), DisplacementT)) -end - -function init_iterates( :: Type{DisplacementT}, dim::ImageSize) where DisplacementT - return init_rest(init_primal(zeros(dim...), DisplacementT)) -end - -############################################## -# Weighting for different displacements types -############################################## - -norm²weight( :: Type{DisplacementConstant}, sz ) = prod(sz) -norm²weight( :: Type{DisplacementFull}, sz ) = 1 - -############ -# Algorithm -############ - -function solve( :: Type{DisplacementT}; - dim :: ImageSize, - iterate = AlgTools.simple_iterate, - params::NamedTuple) where DisplacementT - - ###################### - # Initialise iterates - ###################### - - x, y, Δx, Δy, x̄ = init_iterates(DisplacementT, dim) - init_data = (params.init == :data) - - # … for tracking cumulative movement - if DisplacementT == DisplacementConstant - ucumul = [0.0, 0.0] - else - ucumul = [NaN, NaN] - end - - ############################################# - # Extract parameters and set up step lengths - ############################################# - - α, ρ, λ, θ, T = params.α, params.ρ, params.λ, params.θ, params.timestep - R_K² = max(∇₂_norm₂₂_est², ∇₂_norm₂∞_est²*params.dynrange^2) - γ = min(1, λ*norm²weight(DisplacementT, size(x.x))) - τ, σ, σ̃, ρ̃ = step_lengths(params, γ, R_K²) - - kernel = params.kernel - - #################### - # Run the algorithm - #################### - - b_next_filt=nothing - - v = iterate(params) do verbose :: Function, - b :: Image, - 🚫unused_v_known :: DisplacementT, - b_next :: Image - - #################################### - # Smooth data for Horn–Schunck term - #################################### - - b_filt, b_next_filt = filter_hs(b, b_next, b_next_filt, kernel) - - ############################ - # Construct K for this step - ############################ - - K! = (yʹ, xʹ) -> begin - # Optical flow part - @. yʹ.flow = b_filt - flow!(yʹ.flow, Δx.x, xʹ.u) - @. yʹ.flow = yʹ.flow - b_next_filt - # TV part - ∇₂!(yʹ.tv, xʹ.x) - end - Kᵀ! = (xʹ, yʹ) -> begin - # Optical flow part: ∇b_k ⋅ y - # - # TODO: This really should depend x.u, but x.u is zero. - # - pointwise_gradiprod_2dᵀ!(xʹ.u, yʹ.flow, b_filt) - # TV part - ∇₂ᵀ!(xʹ.x, yʹ.tv) - end - - ################## - # Prediction step - ################## - - if init_data - x .= b - init_data = false - end - - pdflow!(x.x, Δx.x, y.tv, Δy.tv, y.flow, Δy.flow, x.u, params.dual_flow) - - # Predict zero displacement - x.u .= 0 - if params.prox_predict - K!(Δy, x) - @. y.tv = (y.tv + σ̃*Δy.tv)/(1 + σ̃*(ρ̃+ρ/α)) - proj_norm₂₁ball!(y.tv, α) - @. y.flow = (y.flow+σ̃*Δy.flow)/(1+σ̃*(ρ̃+1/θ)) - end - - ############ - # PDPS step - # - # NOTE: For DisplacementConstant, the x.u update is supposed to be with - # respect to the 𝟙^*𝟙 norm/inner product that makes the norm equivalent - # to full-space norm when restricted to constant displacements. Since - # `OpticalFlow.pointwise_gradiprod_2dᵀ!` already uses this inner product, - # and the λ-weighted term in the problem is with respect to this norm, - # all the norm weights disappear in this update. - ############ - - Kᵀ!(Δx, y) # primal step: - @. x̄.x = x.x # | save old x for over-relax - @. x̄.u = x.u # | - @. x.x = (x.x-τ*(Δx.x-b))/(1+τ) # | prox - @. x.u = (x.u-τ*Δx.u)/(1+τ*λ) # | - @. x̄.x = 2x.x - x̄.x # over-relax - @. x̄.u = 2x.u - x̄.u # | - K!(Δy, x̄) # dual step: y - @. y.tv = (y.tv + σ*Δy.tv)/(1 + σ*ρ/α) # | - proj_norm₂₁ball!(y.tv, α) # | prox - @. y.flow = (y.flow+σ*Δy.flow)/(1+σ/θ) - - if DisplacementT == DisplacementConstant - ucumul .+= x.u - end - - ######################################################## - # Give function value and cumulative movement if needed - ######################################################## - v = verbose() do - K!(Δy, x) - value = (norm₂²(b-x.x)/2 + θ*norm₂²(Δy.flow) - + λ*norm₂²(x.u)/2 + α*γnorm₂₁(Δy.tv, ρ)) - - value, x.x, ucumul, nothing - end - - return v - end - - return x, y, v -end - -end # Module - -
--- a/src/AlgorithmDualScaling.jl Thu Apr 25 19:16:17 2024 +0300 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,125 +0,0 @@ -#################################################################### -# 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 ..OpticalFlow: ImageSize, - Image, - pdflow!, - DualScaling - -######################### -# Iterate initialisation -######################### - -function init_rest(x::Image) - imdim=size(x) - - y = zeros(2, imdim...) - Δx = copy(x) - Δy = copy(y) - x̄ = copy(x) - - return x, y, Δx, Δy, x̄ -end - -function init_iterates(xinit::Image) - return init_rest(copy(xinit)) -end - -function init_iterates(dim::ImageSize) - return init_rest(zeros(dim...)) -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.0 - Λ = params.Λ - τ₀, σ₀ = params.τ₀, params.σ₀ - - τ = τ₀/γ - @assert(1+γ*τ ≥ Λ) - σ = σ₀*1/(τ*R_K²) - - println("Step length parameters: τ=$(τ), σ=$(σ)") - - ###################### - # Initialise iterates - ###################### - - x, y, Δx, Δy, x̄ = init_iterates(dim) - init_data = (params.init == :data) - - #################### - # Run the algorithm - #################### - - v = iterate(params) do verbose :: Function, - b :: Image, - v_known :: DisplacementT, - 🚫unused_b_next :: Image - - ################## - # Prediction step - ################## - if init_data - x .= b - init_data = false - end - - pdflow!(x, Δx, y, Δy, v_known, DualScaling()) - - ############ - # PDPS step - ############ - - ∇₂ᵀ!(Δx, y) # primal step: - @. x̄ = x # | save old x for over-relax - @. x = (x-τ*(Δx-b))/(1+τ) # | prox - @. x̄ = 2x - x̄ # over-relax - ∇₂!(Δy, x̄) # dual step: y - @. y = (y + σ*Δy)/(1 + σ*ρ/α) # | - proj_norm₂₁ball!(y, α) # | prox - - ################################ - # Give function value if needed - ################################ - v = verbose() do - ∇₂!(Δy, x) - value = norm₂²(b-x)/2 + params.α*γnorm₂₁(Δy, params.ρ) - value, x, [NaN, NaN], nothing - end - - v - end - - return x, y, v -end - -end # Module - -
--- a/src/AlgorithmGreedy.jl Thu Apr 25 19:16:17 2024 +0300 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,125 +0,0 @@ -#################################################################### -# 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 ..OpticalFlow: ImageSize, - Image, - pdflow!, - Greedy - -######################### -# Iterate initialisation -######################### - -function init_rest(x::Image) - imdim=size(x) - - y = zeros(2, imdim...) - Δx = copy(x) - Δy = copy(y) - x̄ = copy(x) - - return x, y, Δx, Δy, x̄ -end - -function init_iterates(xinit::Image) - return init_rest(copy(xinit)) -end - -function init_iterates(dim::ImageSize) - return init_rest(zeros(dim...)) -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.0 - Λ = params.Λ - τ₀, σ₀ = params.τ₀, params.σ₀ - - τ = τ₀/γ - @assert(1+γ*τ ≥ Λ) - σ = σ₀*1/(τ*R_K²) - - println("Step length parameters: τ=$(τ), σ=$(σ)") - - ###################### - # Initialise iterates - ###################### - - x, y, Δx, Δy, x̄ = init_iterates(dim) - init_data = (params.init == :data) - - #################### - # Run the algorithm - #################### - - v = iterate(params) do verbose :: Function, - b :: Image, - v_known :: DisplacementT, - 🚫unused_b_next :: Image - - ################## - # Prediction step - ################## - if init_data - x .= b - init_data = false - end - - pdflow!(x, Δx, y, Δy, v_known, Greedy()) - - ############ - # PDPS step - ############ - - ∇₂ᵀ!(Δx, y) # primal step: - @. x̄ = x # | save old x for over-relax - @. x = (x-τ*(Δx-b))/(1+τ) # | prox - @. x̄ = 2x - x̄ # over-relax - ∇₂!(Δy, x̄) # dual step: y - @. y = (y + σ*Δy)/(1 + σ*ρ/α) # | - proj_norm₂₁ball!(y, α) # | prox - - ################################ - # Give function value if needed - ################################ - v = verbose() do - ∇₂!(Δy, x) - value = norm₂²(b-x)/2 + params.α*γnorm₂₁(Δy, params.ρ) - 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/AlgorithmNew.jl Thu Apr 25 14:20:38 2024 -0500 @@ -0,0 +1,127 @@ +#################################################################### +# Predictive online PDPS for optical flow with known velocity field +#################################################################### + +__precompile__() + +module AlgorithmNew + +# Default identifier for when none is given by predictor +identifier = "pdps_known_noprediction" + +using Printf + +using AlgTools.Util +import AlgTools.Iterate +using ImageTools.Gradient + +using ..OpticalFlow: ImageSize, + Image, + pdflow! + +######################### +# Iterate initialisation +######################### + +function init_rest(x::Image) + imdim=size(x) + + y = zeros(2, imdim...) + Δx = copy(x) + Δy = copy(y) + x̄ = copy(x) + + return x, y, Δx, Δy, x̄ +end + +function init_iterates(xinit::Image) + return init_rest(copy(xinit)) +end + +function init_iterates(dim::ImageSize) + return init_rest(zeros(dim...)) +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.0 + Λ = params.Λ + τ₀, σ₀ = params.τ₀, params.σ₀ + + τ = τ₀/γ + @assert(1+γ*τ ≥ Λ) + σ = σ₀*1/(τ*R_K²) + + println("Step length parameters: τ=$(τ), σ=$(σ)") + + ###################### + # Initialise iterates + ###################### + + x, y, Δx, Δy, x̄ = init_iterates(dim) + init_data = (params.init == :data) + + #################### + # Run the algorithm + #################### + + v = iterate(params) do verbose :: Function, + b :: Image, + v_known :: DisplacementT, + 🚫unused_b_next :: Image + + ################## + # Prediction step + ################## + if init_data + x .= b + init_data = false + end + + if haskey(params, :predictor) && ~isnothing(params.predictor) + pdflow!(x, Δx, y, Δy, v_known, params.predictor) + end + + ############ + # PDPS step + ############ + + ∇₂ᵀ!(Δx, y) # primal step: + @. x̄ = x # | save old x for over-relax + @. x = (x-τ*(Δx-b))/(1+τ) # | prox + @. x̄ = 2x - x̄ # over-relax + ∇₂!(Δy, x̄) # dual step: y + @. y = (y + σ*Δy)/(1 + σ*ρ/α) # | + proj_norm₂₁ball!(y, α) # | prox + + ################################ + # Give function value if needed + ################################ + v = verbose() do + ∇₂!(Δy, x) + value = norm₂²(b-x)/2 + params.α*γnorm₂₁(Δy, params.ρ) + value, x, [NaN, NaN], nothing + end + + v + end + + return x, y, v +end + +end # Module + +
--- a/src/AlgorithmNoPrediction.jl Thu Apr 25 19:16:17 2024 +0300 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,124 +0,0 @@ -#################################################################### -# 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 ..OpticalFlow: ImageSize, - Image, - pdflow! - -######################### -# Iterate initialisation -######################### - -function init_rest(x::Image) - imdim=size(x) - - y = zeros(2, imdim...) - Δx = copy(x) - Δy = copy(y) - x̄ = copy(x) - - return x, y, Δx, Δy, x̄ -end - -function init_iterates(xinit::Image) - return init_rest(copy(xinit)) -end - -function init_iterates(dim::ImageSize) - return init_rest(zeros(dim...)) -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.0 - Λ = params.Λ - τ₀, σ₀ = params.τ₀, params.σ₀ - - τ = τ₀/γ - @assert(1+γ*τ ≥ Λ) - σ = σ₀*1/(τ*R_K²) - - println("Step length parameters: τ=$(τ), σ=$(σ)") - - ###################### - # Initialise iterates - ###################### - - x, y, Δx, Δy, x̄ = init_iterates(dim) - init_data = (params.init == :data) - - #################### - # Run the algorithm - #################### - - v = iterate(params) do verbose :: Function, - b :: Image, - v_known :: DisplacementT, - 🚫unused_b_next :: Image - - ################## - # Prediction step - ################## - if init_data - x .= b - init_data = false - end - - # No prediction - - ############ - # PDPS step - ############ - - ∇₂ᵀ!(Δx, y) # primal step: - @. x̄ = x # | save old x for over-relax - @. x = (x-τ*(Δx-b))/(1+τ) # | prox - @. x̄ = 2x - x̄ # over-relax - ∇₂!(Δy, x̄) # dual step: y - @. y = (y + σ*Δy)/(1 + σ*ρ/α) # | - proj_norm₂₁ball!(y, α) # | prox - - ################################ - # Give function value if needed - ################################ - v = verbose() do - ∇₂!(Δy, x) - value = norm₂²(b-x)/2 + params.α*γnorm₂₁(Δy, params.ρ) - value, x, [NaN, NaN], nothing - end - - v - end - - return x, y, v -end - -end # Module - -
--- a/src/AlgorithmPrimalOnly.jl Thu Apr 25 19:16:17 2024 +0300 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,124 +0,0 @@ -#################################################################### -# 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 ..OpticalFlow: ImageSize, - Image, - pdflow! - -######################### -# Iterate initialisation -######################### - -function init_rest(x::Image) - imdim=size(x) - - y = zeros(2, imdim...) - Δx = copy(x) - Δy = copy(y) - x̄ = copy(x) - - return x, y, Δx, Δy, x̄ -end - -function init_iterates(xinit::Image) - return init_rest(copy(xinit)) -end - -function init_iterates(dim::ImageSize) - return init_rest(zeros(dim...)) -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.0 - Λ = params.Λ - τ₀, σ₀ = params.τ₀, params.σ₀ - - τ = τ₀/γ - @assert(1+γ*τ ≥ Λ) - σ = σ₀*1/(τ*R_K²) - - println("Step length parameters: τ=$(τ), σ=$(σ)") - - ###################### - # Initialise iterates - ###################### - - x, y, Δx, Δy, x̄ = init_iterates(dim) - init_data = (params.init == :data) - - #################### - # Run the algorithm - #################### - - v = iterate(params) do verbose :: Function, - b :: Image, - v_known :: DisplacementT, - 🚫unused_b_next :: Image - - ################## - # Prediction step - ################## - if init_data - x .= b - init_data = false - end - - pdflow!(x, Δx, y, Δy, v_known, false) - - ############ - # PDPS step - ############ - - ∇₂ᵀ!(Δx, y) # primal step: - @. x̄ = x # | save old x for over-relax - @. x = (x-τ*(Δx-b))/(1+τ) # | prox - @. x̄ = 2x - x̄ # over-relax - ∇₂!(Δy, x̄) # dual step: y - @. y = (y + σ*Δy)/(1 + σ*ρ/α) # | - proj_norm₂₁ball!(y, α) # | prox - - ################################ - # Give function value if needed - ################################ - v = verbose() do - ∇₂!(Δy, x) - value = norm₂²(b-x)/2 + params.α*γnorm₂₁(Δy, params.ρ) - value, x, [NaN, NaN], nothing - end - - v - end - - return x, y, v -end - -end # Module - -
--- a/src/AlgorithmRotation.jl Thu Apr 25 19:16:17 2024 +0300 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,125 +0,0 @@ -#################################################################### -# 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 ..OpticalFlow: ImageSize, - Image, - pdflow!, - Rotation - -######################### -# Iterate initialisation -######################### - -function init_rest(x::Image) - imdim=size(x) - - y = zeros(2, imdim...) - Δx = copy(x) - Δy = copy(y) - x̄ = copy(x) - - return x, y, Δx, Δy, x̄ -end - -function init_iterates(xinit::Image) - return init_rest(copy(xinit)) -end - -function init_iterates(dim::ImageSize) - return init_rest(zeros(dim...)) -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.0 - Λ = params.Λ - τ₀, σ₀ = params.τ₀, params.σ₀ - - τ = τ₀/γ - @assert(1+γ*τ ≥ Λ) - σ = σ₀*1/(τ*R_K²) - - println("Step length parameters: τ=$(τ), σ=$(σ)") - - ###################### - # Initialise iterates - ###################### - - x, y, Δx, Δy, x̄ = init_iterates(dim) - init_data = (params.init == :data) - - #################### - # Run the algorithm - #################### - - v = iterate(params) do verbose :: Function, - b :: Image, - v_known :: DisplacementT, - 🚫unused_b_next :: Image - - ################## - # Prediction step - ################## - if init_data - x .= b - init_data = false - end - - pdflow!(x, Δx, y, Δy, v_known, Rotation()) - - ############ - # PDPS step - ############ - - ∇₂ᵀ!(Δx, y) # primal step: - @. x̄ = x # | save old x for over-relax - @. x = (x-τ*(Δx-b))/(1+τ) # | prox - @. x̄ = 2x - x̄ # over-relax - ∇₂!(Δy, x̄) # dual step: y - @. y = (y + σ*Δy)/(1 + σ*ρ/α) # | - proj_norm₂₁ball!(y, α) # | prox - - ################################ - # Give function value if needed - ################################ - v = verbose() do - ∇₂!(Δy, x) - value = norm₂²(b-x)/2 + params.α*γnorm₂₁(Δy, params.ρ) - value, x, [NaN, NaN], nothing - end - - v - end - - return x, y, v -end - -end # Module - -
--- a/src/AlgorithmZeroDual.jl Thu Apr 25 19:16:17 2024 +0300 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,125 +0,0 @@ -#################################################################### -# Predictive online PDPS for optical flow with known velocity field -#################################################################### - -__precompile__() - -module AlgorithmZeroDual - -identifier = "pdps_known_zerodual" - -using Printf - -using AlgTools.Util -import AlgTools.Iterate -using ImageTools.Gradient - -using ..OpticalFlow: ImageSize, - Image, - pdflow! - -######################### -# Iterate initialisation -######################### - -function init_rest(x::Image) - imdim=size(x) - - y = zeros(2, imdim...) - Δx = copy(x) - Δy = copy(y) - x̄ = copy(x) - - return x, y, Δx, Δy, x̄ -end - -function init_iterates(xinit::Image) - return init_rest(copy(xinit)) -end - -function init_iterates(dim::ImageSize) - return init_rest(zeros(dim...)) -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.0 - Λ = params.Λ - τ₀, σ₀ = params.τ₀, params.σ₀ - - τ = τ₀/γ - @assert(1+γ*τ ≥ Λ) - σ = σ₀*1/(τ*R_K²) - - println("Step length parameters: τ=$(τ), σ=$(σ)") - - ###################### - # Initialise iterates - ###################### - - x, y, Δx, Δy, x̄ = init_iterates(dim) - init_data = (params.init == :data) - - #################### - # Run the algorithm - #################### - - v = iterate(params) do verbose :: Function, - b :: Image, - v_known :: DisplacementT, - 🚫unused_b_next :: Image - - ################## - # Prediction step - ################## - if init_data - x .= b - init_data = false - end - - pdflow!(x, Δx, y, Δy, v_known, false) - y .= zeros(size(y)...) - - ############ - # PDPS step - ############ - - ∇₂ᵀ!(Δx, y) # primal step: - @. x̄ = x # | save old x for over-relax - @. x = (x-τ*(Δx-b))/(1+τ) # | prox - @. x̄ = 2x - x̄ # over-relax - ∇₂!(Δy, x̄) # dual step: y - @. y = (y + σ*Δy)/(1 + σ*ρ/α) # | - proj_norm₂₁ball!(y, α) # | prox - - ################################ - # Give function value if needed - ################################ - v = verbose() do - ∇₂!(Δy, x) - value = norm₂²(b-x)/2 + params.α*γnorm₂₁(Δy, params.ρ) - value, x, [NaN, NaN], nothing - end - - v - end - - return x, y, v -end - -end # Module - -
--- a/src/ImGenerate.jl Thu Apr 25 19:16:17 2024 +0300 +++ b/src/ImGenerate.jl Thu Apr 25 14:20:38 2024 -0500 @@ -100,10 +100,10 @@ θ = 2π*rand(rng,Float64) r = params.shake return [r*cos(θ), r*sin(θ)] - end + end else error("Unknown shaketype $(params.shaketype)") - end + end end pixelwise = (shakefn, sz) -> () -> make_const_u(shakefn(), sz) @@ -181,15 +181,14 @@ :: Type{DisplacementConstant}, datachannel :: Channel{OnlineData{DisplacementConstant}}, params :: NamedTuple) - # Set up counter and zero factor for stabilisation @assert(params.maxiter ≥ maximum(params.stable_interval)) indx = 1 - zero_factor = indx in params.stable_interval ? 0.0 : 1.0 + zero_factor = indx in params.stable_interval ? 0.0 : 1.0 - # Restart the seed to enable comparison across predictors - Random.seed!(rng,951508) # algreadygood + # Restart the seed to enable comparison across predictors + Random.seed!(rng, if haskey(params, :seed) params.seed else 951508 end) nextv = shake(params) v_true = zero_factor.*nextv() @@ -228,18 +227,17 @@ # PETscan ######################################################################## function generate_sinogram(im, sz, - :: Type{DisplacementConstant}, - datachannel :: Channel{PetOnlineData{DisplacementConstant}}, - params :: NamedTuple) + :: Type{DisplacementConstant}, + datachannel :: Channel{PetOnlineData{DisplacementConstant}}, + params :: NamedTuple) # Set up counter and zero factor for stabilisation @assert(params.maxiter ≥ maximum(params.stable_interval)) indx = 1 zero_factor = indx in params.stable_interval ? 0.0 : 1.0 - # Restart the seed to enable comparison across predictors - # Random.seed!(rng,314159) # for shepp - Random.seed!(rng,9182737465) # for brain + # Restart the seed to enable comparison across predictors + Random.seed!(rng, if haskey(params, :seed) params.seed else 314159 end) nextv = shake(params) v_true = zero_factor.*nextv() @@ -283,7 +281,7 @@ # Next theta theta_true = zero_factor*rotatebytheta(params) theta_cumul += theta_true - # Next sinogram mask + # Next sinogram mask S_true = generate_radonmask(params) # Update indx and zero factor indx += 1
--- a/src/OpticalFlow.jl Thu Apr 25 19:16:17 2024 +0300 +++ b/src/OpticalFlow.jl Thu Apr 25 14:20:38 2024 -0500 @@ -42,7 +42,8 @@ HornSchunckData, filter_hs, petpdflow!, - DualScaling, Greedy, Rotation + DualScaling, Greedy, Rotation, ZeroDual, PrimalOnly, + identifier ############################################### # Types (several imported from ImageTools.Translate) @@ -60,9 +61,36 @@ ################################# # Struct for flow ################################# -struct DualScaling end +struct DualScaling + exponent :: Integer + threshold :: Real + DualScaling(e = 10, t = 1e-1) = new(e, t) +end + struct Greedy end struct Rotation end +struct ZeroDual end +struct PrimalOnly end + +function identifier(::DualScaling) + "pdps_known_dualscaling" +end + +function identifier(::Rotation) + "pdps_known_rotation" +end + +function identifier(::Greedy) + "pdps_known_greedy" +end + +function identifier(::ZeroDual) + "pdps_known_zerodual" +end + +function identifier(::PrimalOnly) + "pdps_known_primalonly" +end ################################# # Displacement field based flow @@ -97,7 +125,7 @@ end end -function pdflow!(x, Δx, y, Δy, z, Δz, u, dual_flow; threads=:none) +function pdflow!(x, Δx, y, Δy, z, Δz, u, dual_flow :: Bool; threads=:none) if dual_flow @backgroundif (threads==:outer) begin flow!(x, u, Δx; threads=(threads==:inner)) @@ -114,7 +142,7 @@ # Additional method for Greedy function pdflow!(x, Δx, y, Δy, u, flow :: Greedy; threads=:none) - @assert(size(u)==(2,)) + @assert(size(u)==(2,)) Δx .= x Δy .= y flow!(x, u; threads=(threads==:inner)) @@ -125,11 +153,11 @@ inds = abs.(Dxx) .≤ 1e-1 Dxx[inds] .= 1 DΔx[inds] .= 1 - y .= y.* DΔx ./ Dxx + y .= y.* DΔx ./ Dxx end # Additional method for Rotation -function pdflow!(x, Δx, y, Δy, u, flow :: Rotation; threads=:none) +function pdflow!(x, Δx, y, Δy, u, flow :: Rotation; threads=:none) @assert(size(u)==(2,)) Δx .= x flow!(x, u; threads=(threads==:inner)) @@ -164,17 +192,28 @@ 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) + cm = max(flow.threshold,maximum(cc)) + c = 1 .* (1 .- cc./ cm) .^flow.exponent C[1,:,:] .= c C[2,:,:] .= c - y .= C.*y + y .= C.*y end +function pdflow!(x, Δx, y, Δy, u, :: ZeroDual; threads=:none) + @assert(size(u)==(2,)) + flow!(x, u; threads=(threads==:inner)) + y .= 0.0 +end + +function pdflow!(x, Δx, y, Δy, u, :: PrimalOnly; threads=:none) + @assert(size(u)==(2,)) + flow!(x, u; threads=(threads==:inner)) +end ########################## # PET ########################## + function petflow_interp!(x::AbstractImage, tmp::AbstractImage, u::DisplacementConstant, theta_known::DisplacementConstant; threads = false) tmp .= x @@ -186,7 +225,7 @@ petflow! = petflow_interp! -function petpdflow!(x, Δx, y, Δy, u, theta_known, dual_flow; threads=:none) +function petpdflow!(x, Δx, y, Δy, u, theta_known, dual_flow :: Bool; threads=:none) if dual_flow @backgroundif (threads==:outer) begin petflow!(x, Δx, u, theta_known; threads=(threads==:inner)) @@ -214,7 +253,7 @@ inds = abs.(Dxx) .≤ 1e-2 Dxx[inds] .= 1 DΔx[inds] .= 1 - y .= y.* DΔx ./ Dxx + y .= y.* DΔx ./ Dxx end # Method for dual scaling predictor @@ -243,6 +282,15 @@ end end +function petpdflow!(x, Δx, y, Δy, u, theta_known, :: ZeroDual; threads=:none) + petflow!(x, Δx, u, theta_known, threads=(threads==:inner)) + y .= 0.0 +end + +function petpdflow!(x, Δx, y, Δy, u, theta_known, :: PrimalOnly; threads=:none) + petflow!(x, Δx, u, theta_known, threads=(threads==:inner)) +end + ########################## # Linearised optical flow ########################## @@ -353,7 +401,7 @@ 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. + # instruction usage. g = zeros(2, size(b)...) ∇₂c!(g, b) @inbounds for i=1:size(b, 1) @@ -420,7 +468,7 @@ f = x -> simple_imfilter(x, kernel; threads=true) end - # We already filtered b in the previous step (b_next in that step) + # 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)
--- a/src/OpticalFlow.jl.orig Thu Apr 25 19:16:17 2024 +0300 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,280 +0,0 @@ -################################ -# Code relevant to optical flow -################################ - -__precompile__() - -module OpticalFlow - -using AlgTools.Util -using ImageTools.Gradient -import ImageTools.Translate -using ImageTools.ImFilter - -########## -# 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 - -############################################### -# 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 - -########################## -# 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=false) - 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
--- a/src/PET/AlgorithmDualScaling.jl Thu Apr 25 19:16:17 2024 +0300 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,181 +0,0 @@ -#################################################################### -# 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!, - DualScaling - -######################### -# 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, DualScaling()) - - 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) - println("Step length parameters: L=$(L)") - τ = τ₀/L - σ = σ₀*(1-τ₀)/(R_K²*τ) - end - 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/PET/AlgorithmGreedy.jl Thu Apr 25 19:16:17 2024 +0300 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,181 +0,0 @@ -#################################################################### -# 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!, - Greedy - -######################### -# 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, Greedy()) - - 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) - println("Step length parameters: L=$(L)") - τ = τ₀/L - σ = σ₀*(1-τ₀)/(R_K²*τ) - end - 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/AlgorithmNew.jl Thu Apr 25 14:20:38 2024 -0500 @@ -0,0 +1,183 @@ +#################################################################### +# Predictive online PDPS for optical flow with known velocity field +#################################################################### + +__precompile__() + +module AlgorithmNew + +# Default identifier for when none is given by predictor +identifier = "pdps_known_noprediction" + +using Printf + +using AlgTools.Util +import AlgTools.Iterate +using ImageTools.Gradient +using ImageTools.Translate + +using ImageTransformations +using Images, CoordinateTransformations, Rotations, OffsetArrays +using ImageCore, Interpolations + +using ...Radon +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 + ################### + + if haskey(params, :predictor) && ~isnothing(params.predictor) + petpdflow!(x, Δx, y, Δy, v_known, theta_known, params.predictor) + end + + 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) + println("Step length parameters: L=$(L)") + τ = τ₀/L + σ = σ₀*(1-τ₀)/(R_K²*τ) + end + 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/PET/AlgorithmNoPrediction.jl Thu Apr 25 19:16:17 2024 +0300 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,180 +0,0 @@ -#################################################################### -# 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) - println("Step length parameters: L=$(L)") - τ = τ₀/L - σ = σ₀*(1-τ₀)/(R_K²*τ) - end - 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/PET/AlgorithmPrimalOnly.jl Thu Apr 25 19:16:17 2024 +0300 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,180 +0,0 @@ -#################################################################### -# 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) - println("Step length parameters: L=$(L)") - τ = τ₀/L - σ = σ₀*(1-τ₀)/(R_K²*τ) - end - 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/PET/AlgorithmProximal.jl Thu Apr 25 19:16:17 2024 +0300 +++ b/src/PET/AlgorithmProximal.jl Thu Apr 25 14:20:38 2024 -0500 @@ -15,14 +15,14 @@ using ImageTools.Gradient using ImageTools.Translate -using ..Radon using ImageTransformations using Images, CoordinateTransformations, Rotations, OffsetArrays using ImageCore, Interpolations -using ..OpticalFlow: ImageSize, - Image, - petpdflow! +using ...Radon +using ...OpticalFlow: ImageSize, + Image, + petpdflow! ######################### # Iterate initialisation @@ -103,9 +103,9 @@ iterate = AlgTools.simple_iterate, params::NamedTuple) where DisplacementT - ################################ + ################################ # Extract and set up parameters - ################################ + ################################ α, ρ = params.α, params.ρ R_K² = ∇₂_norm₂₂_est² γ = 1 @@ -139,9 +139,9 @@ v_known :: DisplacementT, theta_known :: DisplacementT, b_true :: Image, - S :: Image + S :: Image - ################### + ################### # Prediction steps ################### @@ -151,14 +151,14 @@ @. oldx = x end - ############################## + ############################## # Proximal step of prediction ############################## ∇₂!(Δy, x) @. y = (y + σ̃*Δy)/(1 + σ̃*(ρ̃+ρ/α)) - #@. cc = y + 1000000*σ̃*Δy + #@. cc = y + 1000000*σ̃*Δy #@. y = (y + σ̃*Δy)/(1 + σ̃*(ρ̃+ρ/α)) + (1 - 1/(1 + ρ̃*σ̃))*cc - proj_norm₂₁ball!(y, α) + proj_norm₂₁ball!(y, α) ############ # PDPS step @@ -177,7 +177,7 @@ ##################### # L update if needed - ##################### + ##################### if params.L_experiment petgrad!(petgradx, x, b, c, S) petgrad!(oldpetgradx, oldx, b, c, S) @@ -186,18 +186,18 @@ println("Step length parameters: L=$(L)") τ = τ₀/L σ = σ₀*(1-τ₀)/(R_K²*τ) - end - end + end + end ################################ # Give function value if needed ################################ - v = verbose() do + v = verbose() do ∇₂!(Δy, x) value = λ*petvalue(x, b, c) + params.α*norm₂₁(Δy) value, x, [NaN, NaN], nothing - end + end v end
--- a/src/PET/AlgorithmRotation.jl Thu Apr 25 19:16:17 2024 +0300 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,181 +0,0 @@ -#################################################################### -# 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!, - Rotation - -######################### -# 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, Rotation()) - - 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) - println("Step length parameters: L=$(L)") - τ = τ₀/L - σ = σ₀*(1-τ₀)/(R_K²*τ) - end - 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/PET/AlgorithmZeroDual.jl Thu Apr 25 19:16:17 2024 +0300 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,181 +0,0 @@ -#################################################################### -# Predictive online PDPS for optical flow with known velocity field -#################################################################### - -__precompile__() - -module AlgorithmZeroDual - -identifier = "pdps_known_zerodual" - -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) - y .= zeros(size(y)...) - - 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) - println("Step length parameters: L=$(L)") - τ = τ₀/L - σ = σ₀*(1-τ₀)/(R_K²*τ) - end - 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/PET/ImGenerate.jl Thu Apr 25 19:16:17 2024 +0300 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,309 +0,0 @@ -################### -# 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(314159) - -# Added for PET -import PoissonRandom: pois_rand -import Random: shuffle -import Images: center, warp, imresize -import CoordinateTransformations: recenter -import Rotations: RotMatrix -import Interpolations: Flat -import MAT: matread - - -############## -# Our exports -############## - -export ImGen, - OnlineData, - imgen_square, - imgen_shake, - PetOnlineData, - imgen_shepplogan_radon, - imgen_brainphantom_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*randn(rng) - return r -end - -function generate_radonmask(params) - imdim = params.radondims - sino_sparse = params.sino_sparsity - numzero = Int64(round(sino_sparse*imdim[1]*imdim[2])) - numone = imdim[1]*imdim[2]-numzero - 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) - - - # Set up counter and zero factor for stabilisation - @assert(params.maxiter ≥ maximum(params.stable_interval)) - indx = 1 - zero_factor = indx in params.stable_interval ? 0.0 : 1.0 - - # Restart the seed to enable comparison across predictors - Random.seed!(rng,951508) # alreadygood - - - nextv = shake(params) - v_true = zero_factor.*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 = zero_factor.*nextv() - v_cumul .+= v_true - # Update indx and zero factor - indx += 1 - zero_factor = indx in params.stable_interval ? 0.0 : 1.0 - 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_sinogram(im, sz, - :: Type{DisplacementConstant}, - datachannel :: Channel{PetOnlineData{DisplacementConstant}}, - params :: NamedTuple) - - # Set up counter and zero factor for stabilisation - @assert(params.maxiter ≥ maximum(params.stable_interval)) - indx = 1 - zero_factor = indx in params.stable_interval ? 0.0 : 1.0 - - # Restart the seed to enable comparison across predictors - Random.seed!(rng,314159) # for shepp - # Random.seed!(rng,9182737465) # for brain - - nextv = shake(params) - v_true = zero_factor.*nextv() - v_cumul = copy(v_true) - - S_true = generate_radonmask(params) - - theta_true = zero_factor*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 = zero_factor.*nextv() - v_cumul .+= v_true - # Next theta - theta_true = zero_factor*rotatebytheta(params) - theta_cumul += theta_true - # Next sinogram mask - S_true = generate_radonmask(params) - # Update indx and zero factor - indx += 1 - zero_factor = indx in params.stable_interval ? 0.0 : 1.0 - end -end - -function imgen_shepplogan_radon(sz) - im = convert(Array{Float64},TestImages.shepp_logan(sz[1], highContrast=true)) - dynrange = maximum(im) - return ImGen(curry(generate_sinogram, im, sz), sz, 1, dynrange, "shepplogan$(sz[1])x$(sz[2])") -end - -function imgen_brainphantom_radon(sz) - data = matread("src/PET/phantom_slice.mat") - im = normalise(imresize(convert(Array{Float64},data["square_data"]),sz)) - dynrange = maximum(im) - return ImGen(curry(generate_sinogram, im, sz), sz, 1, dynrange, "brainphantom$(sz[1])x$(sz[2])") -end - -normalise = (data) -> data./maximum(data) -end # Module
--- a/src/PET/OpticalFlow.jl Thu Apr 25 19:16:17 2024 +0300 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,430 +0,0 @@ -################################ -# 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!, - DualScaling, Greedy, Rotation - -############################################### -# 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} - - -################################# -# Struct for flow -################################# -struct DualScaling end -struct Greedy end -struct Rotation end - -################################# -# 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, flow :: Greedy; 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, Δy, u, flow :: Rotation; 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, Δx, y, Δy, u, flow :: DualScaling; 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, flow :: Greedy; 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 - Dxx = copy(Δy) - DΔx = copy(Δy) - ∇₂!(Dxx, x) - ∇₂!(DΔx, oldx) - inds = abs.(Dxx) .≤ 1e-2 - Dxx[inds] .= 1 - DΔx[inds] .= 1 - y .= y.* DΔx ./ Dxx -end - -# Method for dual scaling predictor -function petpdflow!(x, Δx, y, Δy, u, theta_known, flow :: DualScaling; 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) - cm = max(1e-12,maximum(cc)) - c = 1 .* (1 .- cc./ cm) .^(10) - C[1,:,:] .= c - C[2,:,:] .= c - y .= C.*y -end - -# Method for rotation prediction (exploiting property of inverse rotation) -function petpdflow!(x, Δx, y, Δy, u, theta_known, flow :: Rotation; threads=:none) - @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 -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
--- a/src/PET/PET.jl Thu Apr 25 19:16:17 2024 +0300 +++ b/src/PET/PET.jl Thu Apr 25 14:20:38 2024 -0500 @@ -18,7 +18,7 @@ using AlgTools.StructTools using AlgTools.LinkedLists using AlgTools.Comms -using ImageTools.Visualise: secs_ns, grayimg, do_visualise +using ImageTools.Visualise: secs_ns, grayimg, do_visualise using ImageTools.ImFilter: gaussian # For PET @@ -26,66 +26,36 @@ ##################### # Load local modules -##################### -include("OpticalFlow.jl") -include("Radon.jl") -include("ImGenerate.jl") -include("AlgorithmDualScaling.jl") -include("AlgorithmGreedy.jl") -include("AlgorithmNoPrediction.jl") -include("AlgorithmPrimalOnly.jl") +#####################a +include("AlgorithmNew.jl") include("AlgorithmProximal.jl") -include("AlgorithmRotation.jl") -include("AlgorithmZeroDual.jl") -include("PlotResults.jl") +#include("PlotResults.jl") -import .AlgorithmDualScaling -import .AlgorithmGreedy -import .AlgorithmNoPrediction -import .AlgorithmPrimalOnly +import .AlgorithmNew import .AlgorithmProximal -import .AlgorithmRotation -import .AlgorithmZeroDual -using .Radon: backproject! -using .ImGenerate -using .OpticalFlow: DisplacementFull, DisplacementConstant -using .PlotResults +using ..Radon: backproject! +using ..ImGenerate +using ..OpticalFlow +using ..Run +#using .PlotResults ############## # Our exports ############## -export demo_petS1, demo_petS2, demo_petS3, +export demo_petS1, demo_petS2, demo_petS3, demo_petS4, demo_petS5, demo_petS6, demo_petS7, - demo_petB1, demo_petB2, demo_petB3, + demo_petB1, demo_petB2, demo_petB3, demo_petB4, demo_petB5, demo_petB6, demo_petB7, - batchrun_shepplogan, batchrun_brainphantom, batchrun_pet, - plot_pet + batchrun_shepplogan, batchrun_brainphantom, batchrun_pet + #plot_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 = ( @@ -107,24 +77,24 @@ stable_interval = Set(0), ) -const p_known₀_pet = ( +const p_known₀_pet = default_params ⬿ ( noise_level = 0.5, shake_noise_level = 0.25, - shake = 1.0, - rotation_factor = 0.15, - rotation_noise_level = 0.035, + shake = 1.0, + rotation_factor = 0.15, + rotation_noise_level = 0.035, α = 0.15, ρ̃₀ = 1.0, σ̃₀ = 1.0, δ = 0.9, σ₀ = 1.0, τ₀ = 0.9, - λ = 1, + λ = 1, radondims = [128,64], sz = (256,256), scale = 1, c = 1.0, - sino_sparsity = 0.5, + sino_sparsity = 0.5, L = 300.0, L_experiment = false, #stable_interval = Set(0), @@ -132,41 +102,43 @@ ) const shepplogan = imgen_shepplogan_radon(p_known₀_pet.sz) +const p_known₀_pets = p_known₀_pet ⬿ ( seed = 314159, ) +const p_known₀_petb = p_known₀_pet ⬿ ( seed = 9182737465, ) const brainphantom = imgen_brainphantom_radon(p_known₀_pet.sz) const shepplogan_experiments_pdps_known = ( - Experiment(AlgorithmDualScaling, DisplacementConstant, shepplogan, - p_known₀_pet), - Experiment(AlgorithmGreedy, DisplacementConstant, shepplogan, - p_known₀_pet), - Experiment(AlgorithmNoPrediction, DisplacementConstant, shepplogan, - p_known₀_pet), - Experiment(AlgorithmPrimalOnly, DisplacementConstant, shepplogan, - p_known₀_pet), + Experiment(AlgorithmNew, DisplacementConstant, shepplogan, + p_known₀_pets ⬿ (predictor=DualScaling(),)), + Experiment(AlgorithmNew, DisplacementConstant, shepplogan, + p_known₀_pets ⬿ (predictor=Greedy(),)), + Experiment(AlgorithmNew, DisplacementConstant, shepplogan, + p_known₀_pets ⬿ (predictor=nothing,),), + Experiment(AlgorithmNew, DisplacementConstant, shepplogan, + p_known₀_pets ⬿ (predictor=PrimalOnly(),)), Experiment(AlgorithmProximal, DisplacementConstant, shepplogan, - p_known₀_pet ⬿ (phantom_ρ = 100,)), - Experiment(AlgorithmRotation, DisplacementConstant, shepplogan, - p_known₀_pet), - Experiment(AlgorithmZeroDual, DisplacementConstant, shepplogan, - p_known₀_pet), + p_known₀_pets ⬿ (phantom_ρ = 100,)), + Experiment(AlgorithmNew, DisplacementConstant, shepplogan, + p_known₀_pets ⬿ (predictor=Rotation(),)), + Experiment(AlgorithmNew, DisplacementConstant, shepplogan, + p_known₀_pets ⬿ (predictor=ZeroDual(),)), ) const brainphantom_experiments_pdps_known = ( - Experiment(AlgorithmDualScaling, DisplacementConstant, brainphantom, - p_known₀_pet), - Experiment(AlgorithmGreedy, DisplacementConstant, brainphantom, - p_known₀_pet), - Experiment(AlgorithmNoPrediction, DisplacementConstant, brainphantom, - p_known₀_pet), - Experiment(AlgorithmPrimalOnly, DisplacementConstant, brainphantom, - p_known₀_pet), + Experiment(AlgorithmNew, DisplacementConstant, brainphantom, + p_known₀_petb ⬿ (predictor=DualScaling(),)), + Experiment(AlgorithmNew, DisplacementConstant, brainphantom, + p_known₀_petb ⬿ (predictor=Greedy(),)), + Experiment(AlgorithmNew, DisplacementConstant, brainphantom, + p_known₀_petb ⬿ (predictor=nothing,),), + Experiment(AlgorithmNew, DisplacementConstant, brainphantom, + p_known₀_petb ⬿ (predictor=PrimalOnly(),)), Experiment(AlgorithmProximal, DisplacementConstant, brainphantom, - p_known₀_pet ⬿ (phantom_ρ = 100,)), - Experiment(AlgorithmRotation, DisplacementConstant, brainphantom, - p_known₀_pet), - Experiment(AlgorithmZeroDual, DisplacementConstant, brainphantom, - p_known₀_pet), + p_known₀_petb ⬿ (phantom_ρ = 100,)), + Experiment(AlgorithmNew, DisplacementConstant, brainphantom, + p_known₀_petb ⬿ (predictor=Rotation(),)), + Experiment(AlgorithmNew, DisplacementConstant, brainphantom, + p_known₀_petb ⬿ (predictor=ZeroDual(),)), ) @@ -178,167 +150,15 @@ brainphantom_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(100*p.α))_$(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_prefix=default_save_prefix, + visfn=iterate_visualise_pet, + datatype=PetOnlineData, save_results=false, save_images=false, visualise=true, @@ -367,6 +187,9 @@ function batchrun_shepplogan(;kwargs...) run_experiments(;experiments=shepplogan_experiments_all, + visfn=iterate_visualise_pet, + datatype=PetOnlineData, + save_prefix=default_save_prefix, save_results=true, save_images=true, visualise=false, @@ -376,6 +199,9 @@ function batchrun_brainphantom(;kwargs...) run_experiments(;experiments=brainphantom_experiments_all, + visfn=iterate_visualise_pet, + datatype=PetOnlineData, + save_prefix=default_save_prefix, save_results=true, save_images=true, visualise=false, @@ -395,21 +221,21 @@ 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) + 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 +function iterate_visualise_pet(datachannel::Channel{PetOnlineData{DisplacementT}}, + st :: State, + step :: Function, + params :: NamedTuple) where DisplacementT try sc = nothing d = take!(datachannel) - for iter=1:params.maxiter + 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 @@ -435,7 +261,7 @@ entry = LogEntry(iter, tm, value, #sc*d.v_cumul_true[1], #sc*d.v_cumul_true[2], - #sc*v[1], sc*v[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), @@ -509,39 +335,6 @@ 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) @@ -557,22 +350,13 @@ # Plotting SSIM and PSNR ######################### -function plot_pet(kwargs...) - ssim_plot("shepplogan") - psnr_plot("shepplogan") - fv_plot("shepplogan") - ssim_plot("brainphantom") - psnr_plot("brainphantom") - fv_plot("brainphantom") -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}}}) +#function plot_pet(kwargs...) +# ssim_plot("shepplogan") +# psnr_plot("shepplogan") +# fv_plot("shepplogan") +# ssim_plot("brainphantom") +# psnr_plot("brainphantom") +# fv_plot("brainphantom") +#end end # Module
--- a/src/PET/Radon.jl Thu Apr 25 19:16:17 2024 +0300 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,85 +0,0 @@ -# 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/PlotResults.jl Thu Apr 25 19:16:17 2024 +0300 +++ b/src/PlotResults.jl Thu Apr 25 14:20:38 2024 -0500 @@ -1,5 +1,3 @@ -__precompile__() - module PlotResults @@ -7,13 +5,12 @@ # Load external modules ######################## -using DelimitedFiles, CSV, DataFrames +using CSV, DataFrames using PlotlyJS using Colors -using XLSX: writetable using Statistics -export fv_plot, ssim_plot, psnr_plot, calculate_statistics +export fv_plot, ssim_plot, psnr_plot global mystart = 38 global myend = 135 @@ -656,66 +653,4 @@ end -function calculate_statistics() - ImName = ("lighthouse200x300", "shepplogan256x256", "brainphantom256x256") - AlgName = ("dualscaling", "greedy","noprediction","primalonly","proximal","rotation","zerodual") - mystart = 41 # Corresponds to the 500th iterate - - # Define an array to store results - results = DataFrame(experiment = String[], α = Float64[], algorithm = String[], psnr_mean1 = Float64[], psnr_mean500 = Float64[], psnr_ci = String[], ssim_mean1 = Float64[], ssim_mean500 = Float64[], ssim_ci = String[]) - - for imname in ImName - for algname in AlgName - directory_path = "./img/" - files = readdir(directory_path) - filtered_files = filter(file -> startswith(file, "$(imname)_pdps_known_$(algname)") && endswith(file, ".txt"), files) - - for file in filtered_files - filename = directory_path * file - data = CSV.File(filename, delim='\t', header=2) |> DataFrame - - # Extract α from filename - α, _, _ = extract_parameters(filename) - - - # Extract SSIM and PSNR columns starting from 1st iteration - ssim_values1 = Float64.(data[:, :ssim]) - psnr_values1 = Float64.(data[:, :psnr]) - - # Extract SSIM and PSNR columns starting from 500th iteration - ssim_values500 = Float64.(data[mystart:end, :ssim]) - psnr_values500 = Float64.(data[mystart:end, :psnr]) - - # Calculate mean and confidence intervals - ssim_mean1 = round(mean(ssim_values1), digits=4) - psnr_mean1 = round(mean(psnr_values1), digits=4) - - ssim_mean500 = round(mean(ssim_values500), digits=4) - psnr_mean500 = round(mean(psnr_values500), digits=4) - ssim_std500 = round(std(ssim_values500), digits=4) - psnr_std500 = round(std(psnr_values500), digits=4) - n = length(ssim_values500) - - ssim_ci_lower = round(ssim_mean500 - 1.96 * ssim_std500 / sqrt(n), digits=4) - ssim_ci_upper = round(ssim_mean500 + 1.96 * ssim_std500 / sqrt(n), digits=4) - psnr_ci_lower = round(psnr_mean500 - 1.96 * psnr_std500 / sqrt(n), digits=4) - psnr_ci_upper = round(psnr_mean500 + 1.96 * psnr_std500 / sqrt(n), digits=4) - - ssim_ci = "$(ssim_ci_lower) - $(ssim_ci_upper)" - psnr_ci = "$(psnr_ci_lower) - $(psnr_ci_upper)" - experiment = "$(imname)" - algorithm = "$(algname)" - - # Append results to DataFrame - push!(results, (experiment, α, algorithm, psnr_mean1, psnr_mean500, psnr_ci, ssim_mean1, ssim_mean500, ssim_ci)) - end - end - end - sort!(results, [:experiment, :α]) - excel_path = "./img/summarystats.xlsx" - if isfile(excel_path) - rm(excel_path) - end - writetable(excel_path, sheetname="Experiments", results) -end end # Module \ No newline at end of file
--- a/src/PredictPDPS.jl Thu Apr 25 19:16:17 2024 +0300 +++ b/src/PredictPDPS.jl Thu Apr 25 14:20:38 2024 -0500 @@ -10,20 +10,8 @@ # Load external modules ######################## -using Printf -using FileIO -#using JLD2 -using Setfield -using ImageQualityIndexes: assess_psnr, assess_ssim -using DelimitedFiles -import GR - +using ImageTools.ImFilter: gaussian using AlgTools.Util -using AlgTools.StructTools -using AlgTools.LinkedLists -using AlgTools.Comms -using ImageTools.Visualise: secs_ns, grayimg, do_visualise -using ImageTools.ImFilter: gaussian ##################### # Load local modules @@ -32,48 +20,29 @@ include("OpticalFlow.jl") include("Radon.jl") include("ImGenerate.jl") -include("Algorithm.jl") -include("AlgorithmBoth.jl") -include("AlgorithmBothGreedyV.jl") -include("AlgorithmBothCumul.jl") +include("Run.jl") +include("AlgorithmProximal.jl") include("AlgorithmBothMulti.jl") -include("AlgorithmBothNL.jl") include("AlgorithmFB.jl") include("AlgorithmFBDual.jl") -include("PlotResults.jl") - - -# Additional -include("AlgorithmProximal.jl") -include("AlgorithmGreedy.jl") -include("AlgorithmRotation.jl") -include("AlgorithmNoPrediction.jl") -include("AlgorithmPrimalOnly.jl") -include("AlgorithmDualScaling.jl") -include("AlgorithmZeroDual.jl") +include("AlgorithmNew.jl") +include("Stats.jl") +#include("PlotResults.jl") include("PET/PET.jl") -import .Algorithm, - .AlgorithmBoth, - .AlgorithmBothGreedyV, - .AlgorithmBothCumul, - .AlgorithmBothMulti, - .AlgorithmBothNL, +import .AlgorithmBothMulti, .AlgorithmFB, .AlgorithmFBDual, .AlgorithmProximal, - .AlgorithmGreedy, - .AlgorithmRotation, - .AlgorithmNoPrediction, - .AlgorithmPrimalOnly, - .AlgorithmDualScaling, - .AlgorithmZeroDual + .AlgorithmNew using .ImGenerate -using .OpticalFlow: DisplacementFull, DisplacementConstant -using .PlotResults +using .OpticalFlow +using .Stats +#using .PlotResults using .PET +using .Run ############## # Our exports @@ -85,37 +54,20 @@ demo_unknown1,demo_unknown2,demo_unknown3, batchrun_denoising, batchrun_predictors, - demo_denoising1, demo_denoising2, demo_denoising3, + demo_denoising1, demo_denoising2, demo_denoising3, demo_denoising4, demo_denoising5, demo_denoising6, demo_denoising7, - demo_petS1, demo_petS2, demo_petS3, + demo_petS1, demo_petS2, demo_petS3, demo_petS4, demo_petS5, demo_petS6, demo_petS7, - demo_petB1, demo_petB2, demo_petB3, + demo_petB1, demo_petB2, demo_petB3, demo_petB4, demo_petB5, demo_petB6, demo_petB7, batchrun_shepplogan, batchrun_brainphantom, batchrun_pet, - plot_denoising, plot_pet, calculate_statistics + calculate_statistics + #plot_denoising, plot_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 = ( @@ -130,8 +82,7 @@ 1000, 2000, 2500, 3000, 4000, 5000, 6000, 7000, 7500, 8000, 9000, 10000, 8700]), pixelwise_displacement=false, - dual_flow = true, - prox_predict = true, + dual_flow = true, # For AlgorithmProximalfrom 2019 paper handle_interrupt = true, init = :zero, plot_movement = false, @@ -141,7 +92,7 @@ const square = imgen_square((200, 300)) const lighthouse = imgen_shake("lighthouse", (200, 300)) -const p_known₀_denoising = ( +const p_known₀_denoising = default_params ⬿ ( noise_level = 0.5, shake_noise_level = 0.025, shake = 2.0, @@ -155,7 +106,7 @@ stable_interval = union(Set(2500:5000),Set(8700:10000)), ) -const p_known₀ = ( +const p_known₀ = default_params ⬿ ( noise_level = 0.5, shake_noise_level = 0.05, shake = 2, @@ -167,7 +118,7 @@ τ₀ = 0.01, ) -const p_unknown₀ = ( +const p_unknown₀ = default_params ⬿ ( noise_level = 0.3, shake_noise_level = 0.05, shake = 2, @@ -185,12 +136,14 @@ ) +# Experiments for 2019 paper + const experiments_pdps_known = ( - Experiment(Algorithm, DisplacementConstant, lighthouse, + Experiment(AlgorithmProximal, DisplacementConstant, lighthouse, p_known₀_denoising ⬿ (phantom_ρ = 0,)), - Experiment(Algorithm, DisplacementConstant, lighthouse, + Experiment(AlgorithmProximal, DisplacementConstant, lighthouse, p_known₀ ⬿ (phantom_ρ = 100,)), - Experiment(Algorithm, DisplacementConstant, square, + Experiment(AlgorithmProximal, DisplacementConstant, square, p_known₀ ⬿ (phantom_ρ = 0,)) ) @@ -214,182 +167,29 @@ experiments_fb_known )) +# Image stabilisation experiments for 2024 paper. PET experiments are in PET/PET.jl + 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(AlgorithmNew, DisplacementConstant, lighthouse, + p_known₀_denoising ⬿ (predictor=DualScaling(),)), + Experiment(AlgorithmNew, DisplacementConstant, lighthouse, + p_known₀_denoising ⬿ (predictor=Greedy(),)), + Experiment(AlgorithmNew, DisplacementConstant, lighthouse, + p_known₀_denoising ⬿ (predictor=nothing,),), + Experiment(AlgorithmNew, DisplacementConstant, lighthouse, + p_known₀_denoising ⬿ (predictor=PrimalOnly(),)), Experiment(AlgorithmProximal, DisplacementConstant, lighthouse, p_known₀_denoising ⬿ (phantom_ρ = 100,)), - Experiment(AlgorithmRotation, DisplacementConstant, lighthouse, - p_known₀_denoising), - Experiment(AlgorithmZeroDual, DisplacementConstant, lighthouse, - p_known₀_denoising), + Experiment(AlgorithmNew, DisplacementConstant, lighthouse, + p_known₀_denoising ⬿ (predictor=Rotation(),)), + Experiment(AlgorithmNew, DisplacementConstant, lighthouse, + p_known₀_denoising ⬿ (predictor=ZeroDual(),)), ) const denoising_experiments_all = Iterators.flatten(( denoising_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(100*p.α))_$(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{OnlineData{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 ####################### @@ -398,6 +198,7 @@ run_experiments(;experiments=(experiment,), save_results=false, save_images=false, + save_prefix=default_save_prefix, visualise=true, recalculate=true, verbose_iter=50, @@ -423,6 +224,7 @@ function batchrun_article(kwargs...) run_experiments(;experiments=experiments_all, + save_prefix=default_save_prefix, save_results=true, save_images=true, visualise=false, @@ -432,6 +234,7 @@ function batchrun_denoising(;kwargs...) run_experiments(;experiments=denoising_experiments_all, + save_prefix=default_save_prefix, save_results=true, save_images=true, visualise=false, @@ -445,167 +248,14 @@ batchrun_pet(;kwargs...) end -###################################################### -# Iterator that does visualisation and log collection -###################################################### - -function iterate_visualise(datachannel::Channel{OnlineData{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.b_noisy, d.v, dnext.b_noisy) 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, ((d.b_noisy, 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)" - # save(File(format"PNG", fn("true", "png")), grayimg(d.b_true)) - # save(File(format"PNG", fn("data", "png")), grayimg(d.b_noisy)) - save(File(format"PNG", fn("reco", "png")), grayimg(x)) - 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 - ######################### # Plotting SSIM and PSNR ######################### -function plot_denoising(kwargs...) - ssim_plot("lighthouse") - psnr_plot("lighthouse") - fv_plot("lighthouse") -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}}}) +#function plot_denoising(kwargs...) +# ssim_plot("lighthouse") +# psnr_plot("lighthouse") +# fv_plot("lighthouse") +#end end # Module
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/src/Run.jl Thu Apr 25 14:20:38 2024 -0500 @@ -0,0 +1,360 @@ +################## +# Experiment running and interactive visualisation +################## + +module Run + +import GR +using Setfield +using Printf +using ImageQualityIndexes: assess_psnr, assess_ssim +using FileIO +#using JLD2 +using DelimitedFiles + +using AlgTools.Util +using AlgTools.StructTools +using AlgTools.LinkedLists +using AlgTools.Comms + +using ImageTools.Visualise: secs_ns, grayimg, do_visualise + +using ..ImGenerate +using ..OpticalFlow: identifier, DisplacementFull, DisplacementConstant + +export run_experiments, + Experiment, + State, + LogEntry, + LogEntryHiFi + +################ +# Experiment +################ + +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 + +################ +# 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 + id = if haskey(p, :predictor) && ~isnothing(p.predictor) + identifier(p.predictor) + else + e.mod.identifier + end + # return "$(ig.name)_$(e.mod.identifier)_$(@sprintf "%x" hash(p))" + return "$(ig.name)_$(id)_$(Int64(100*p.α))_$(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, + visfn = iterate_visualise, + datatype = OnlineData, + recalculate=true, + experiments, + 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 = 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{datatype{e.DisplacementT}}(2) + gentask = Threads.@spawn e.imgen.f(e.DisplacementT, datachannel, e_params) + bind(datachannel, gentask) + + # Run algorithm + iterate = curry(visfn, 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 + + +###################################################### +# Iterator that does visualisation and log collection +###################################################### + +function iterate_visualise(datachannel::Channel{OnlineData{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.b_noisy, d.v, dnext.b_noisy) 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, ((d.b_noisy, 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)" + # save(File(format"PNG", fn("true", "png")), grayimg(d.b_true)) + # save(File(format"PNG", fn("data", "png")), grayimg(d.b_noisy)) + save(File(format"PNG", fn("reco", "png")), grayimg(x)) + 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 + +end # module
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/src/Stats.jl Thu Apr 25 14:20:38 2024 -0500 @@ -0,0 +1,78 @@ +__precompile__() + +module Stats + + +######################## +# Load external modules +######################## + +using CSV, DataFrames +using Statistics + +export calculate_statistics + +function calculate_statistics() + ImName = ("lighthouse200x300", "shepplogan256x256", "brainphantom256x256") + AlgName = ("dualscaling", "greedy","noprediction","primalonly","proximal","rotation","zerodual") + mystart = 41 # Corresponds to the 500th iterate + + # Define an array to store results + results = DataFrame(experiment = String[], α = Float64[], algorithm = String[], psnr_mean1 = Float64[], psnr_mean500 = Float64[], psnr_ci = String[], ssim_mean1 = Float64[], ssim_mean500 = Float64[], ssim_ci = String[]) + + for imname in ImName + for algname in AlgName + directory_path = "./img/" + files = readdir(directory_path) + filtered_files = filter(file -> startswith(file, "$(imname)_pdps_known_$(algname)") && endswith(file, ".txt"), files) + + for file in filtered_files + filename = directory_path * file + data = CSV.File(filename, delim='\t', header=2) |> DataFrame + + # Extract α from filename + α, _, _ = extract_parameters(filename) + + + # Extract SSIM and PSNR columns starting from 1st iteration + ssim_values1 = Float64.(data[:, :ssim]) + psnr_values1 = Float64.(data[:, :psnr]) + + # Extract SSIM and PSNR columns starting from 500th iteration + ssim_values500 = Float64.(data[mystart:end, :ssim]) + psnr_values500 = Float64.(data[mystart:end, :psnr]) + + # Calculate mean and confidence intervals + ssim_mean1 = round(mean(ssim_values1), digits=4) + psnr_mean1 = round(mean(psnr_values1), digits=4) + + ssim_mean500 = round(mean(ssim_values500), digits=4) + psnr_mean500 = round(mean(psnr_values500), digits=4) + ssim_std500 = round(std(ssim_values500), digits=4) + psnr_std500 = round(std(psnr_values500), digits=4) + n = length(ssim_values500) + + ssim_ci_lower = round(ssim_mean500 - 1.96 * ssim_std500 / sqrt(n), digits=4) + ssim_ci_upper = round(ssim_mean500 + 1.96 * ssim_std500 / sqrt(n), digits=4) + psnr_ci_lower = round(psnr_mean500 - 1.96 * psnr_std500 / sqrt(n), digits=4) + psnr_ci_upper = round(psnr_mean500 + 1.96 * psnr_std500 / sqrt(n), digits=4) + + ssim_ci = "$(ssim_ci_lower) - $(ssim_ci_upper)" + psnr_ci = "$(psnr_ci_lower) - $(psnr_ci_upper)" + experiment = "$(imname)" + algorithm = "$(algname)" + + # Append results to DataFrame + push!(results, (experiment, α, algorithm, psnr_mean1, psnr_mean500, psnr_ci, ssim_mean1, ssim_mean500, ssim_ci)) + end + end + end + sort!(results, [:experiment, :α]) + csv_path = "./img/summarystats.csv" + if isfile(csv_path) + rm(csv_path) + end + CSV.write(csv_path, results) +end + +end
