Mercurial > repos > PredictPDPS / changeset
changeset
Initialise independent repo
Tue, 07 Apr 2020 14:19:48 -0500
- author
- Tuomo Valkonen <tuomov@iki.fi>
- date
- Tue, 07 Apr 2020 14:19:48 -0500
- changeset 0
- a55e35d20336
- child 1
- dc846d506f2b
Initialise independent repo
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/Manifest.toml Tue Apr 07 14:19:48 2020 -0500 @@ -0,0 +1,391 @@ +# This file is machine-generated - editing it directly is not advised + +[[AbstractFFTs]] +deps = ["LinearAlgebra"] +git-tree-sha1 = "051c95d6836228d120f5f4b984dd5aba1624f716" +uuid = "621f4979-c628-5d54-868e-fcf4e3e8185c" +version = "0.5.0" + +[[AlgTools]] +deps = ["DelimitedFiles", "Printf"] +path = "../AlgTools" +uuid = "c46e2e78-5339-41fd-a966-983ff60ab8e7" +version = "0.1.0" + +[[AxisArrays]] +deps = ["Dates", "IntervalSets", "IterTools", "RangeArrays"] +git-tree-sha1 = "d63ba0315a1d287c9467e61e932578f2fdd048e0" +uuid = "39de3d68-74b9-583c-8d2d-e117c070f3a9" +version = "0.3.3" + +[[Base64]] +uuid = "2a0f44e3-6c83-55bd-87e4-b1978d98bd5f" + +[[BinaryProvider]] +deps = ["Libdl", "SHA"] +git-tree-sha1 = "5b08ed6036d9d3f0ee6369410b830f8873d4024c" +uuid = "b99e7846-7c00-51b0-8f62-c81ae34c0232" +version = "0.5.8" + +[[CatIndices]] +deps = ["CustomUnitRanges", "OffsetArrays", "Test"] +git-tree-sha1 = "254cf73ea369d2e39bfd6c5eb27a2296cfaed68c" +uuid = "aafaddc9-749c-510e-ac4f-586e18779b91" +version = "0.2.0" + +[[ColorTypes]] +deps = ["FixedPointNumbers", "Random"] +git-tree-sha1 = "7b62b728a5f3dd6ee3b23910303ccf27e82fad5e" +uuid = "3da002f7-5984-5a60-b8a6-cbb66c0b333f" +version = "0.8.1" + +[[ColorVectorSpace]] +deps = ["ColorTypes", "Colors", "FixedPointNumbers", "LinearAlgebra", "SpecialFunctions", "Statistics", "StatsBase"] +git-tree-sha1 = "59a561934dc2ad6aba4f500aa839c0f1b8fa8314" +uuid = "c3611d14-8923-5661-9e6a-0046d554d3a4" +version = "0.8.2" + +[[Colors]] +deps = ["ColorTypes", "FixedPointNumbers", "InteractiveUtils", "Printf", "Reexport"] +git-tree-sha1 = "5ae09fba6989add2f55282453b46d68398f1c5ec" +uuid = "5ae59095-9a9b-59fe-a467-6f913c188581" +version = "0.10.0" + +[[Compat]] +deps = ["Base64", "Dates", "DelimitedFiles", "Distributed", "InteractiveUtils", "LibGit2", "Libdl", "LinearAlgebra", "Markdown", "Mmap", "Pkg", "Printf", "REPL", "Random", "Serialization", "SharedArrays", "Sockets", "SparseArrays", "Statistics", "Test", "UUIDs", "Unicode"] +git-tree-sha1 = "3819f476b6b37ef8ea837070ed831b4ebadfa1e9" +uuid = "34da2185-b29b-5c13-b0c7-acf172513d20" +version = "3.2.0" + +[[ComputationalResources]] +deps = ["Test"] +git-tree-sha1 = "89e7e7ed20af73d9f78877d2b8d1194e7b6ff13d" +uuid = "ed09eef8-17a6-5b46-8889-db040fac31e3" +version = "0.3.0" + +[[ConstructionBase]] +git-tree-sha1 = "a2a6a5fea4d6f730ec4c18a76d27ec10e8ec1c50" +uuid = "187b0558-2788-49d3-abe0-74a17ed4e7c9" +version = "1.0.0" + +[[CustomUnitRanges]] +deps = ["Test"] +git-tree-sha1 = "0a106457a1831555857e18ac9617279c22fc393b" +uuid = "dc8bdbbb-1ca9-579f-8c36-e416f6a65cce" +version = "0.2.0" + +[[DataAPI]] +git-tree-sha1 = "674b67f344687a88310213ddfa8a2b3c76cc4252" +uuid = "9a962f9c-6df0-11e9-0e5d-c546b8b5ee8a" +version = "1.1.0" + +[[DataStructures]] +deps = ["InteractiveUtils", "OrderedCollections"] +git-tree-sha1 = "a1b652fb77ae8ca7ea328fa7ba5aa151036e5c10" +uuid = "864edb3b-99cc-5e75-8d2d-829cb0a9cfe8" +version = "0.17.6" + +[[Dates]] +deps = ["Printf"] +uuid = "ade2ca70-3891-5945-98fb-dc099432e06a" + +[[DelimitedFiles]] +deps = ["Mmap"] +uuid = "8bb1440f-4735-579b-a4ab-409b98df4dab" + +[[Distances]] +deps = ["LinearAlgebra", "Statistics"] +git-tree-sha1 = "23717536c81b63e250f682b0e0933769eecd1411" +uuid = "b4f34e82-e78d-54a5-968a-f98e89d6e8f7" +version = "0.8.2" + +[[Distributed]] +deps = ["Random", "Serialization", "Sockets"] +uuid = "8ba89e20-285c-5b6f-9357-94700520ee1b" + +[[FFTViews]] +deps = ["CustomUnitRanges", "FFTW"] +git-tree-sha1 = "1ba17a68391b2b3a0626ffb597c45e482cf7b5f0" +uuid = "4f61f5a4-77b1-5117-aa51-3ab5ef4ef0cd" +version = "0.3.0" + +[[FFTW]] +deps = ["AbstractFFTs", "FFTW_jll", "IntelOpenMP_jll", "Libdl", "LinearAlgebra", "MKL_jll", "Reexport"] +git-tree-sha1 = "109d82fa4b00429f9afcce873e9f746f11f018d3" +uuid = "7a1cc6ca-52ef-59f5-83cd-3a7055c09341" +version = "1.2.0" + +[[FFTW_jll]] +deps = ["Libdl", "Pkg"] +git-tree-sha1 = "05674f209a6e3387dd103a945b0113eeb64b1a58" +uuid = "f5851436-0d7a-5f13-b9de-f02708fd171a" +version = "3.3.9+3" + +[[FileIO]] +deps = ["Pkg"] +git-tree-sha1 = "80c17c711c41416eb0ac68347dc036be68b37682" +uuid = "5789e2e9-d7fb-5bc7-8068-2c6fae9b9549" +version = "1.2.0" + +[[FixedPointNumbers]] +git-tree-sha1 = "bd1386f890e172ef38e1c735cda58cbf004a7c9a" +uuid = "53c48c17-4a7d-5ca2-90c5-79b7896eea93" +version = "0.7.0" + +[[Future]] +deps = ["Random"] +uuid = "9fa8497b-333b-5362-9e8d-4d0656e87820" + +[[GR]] +deps = ["Base64", "DelimitedFiles", "LinearAlgebra", "Printf", "Random", "Serialization", "Sockets", "Test"] +git-tree-sha1 = "c690c2ab22ac9ee323d9966deae61a089362b25c" +uuid = "28b8d3ca-fb5f-59d9-8090-bfdbd6d07a71" +version = "0.44.0" + +[[Graphics]] +deps = ["Colors", "LinearAlgebra", "NaNMath"] +git-tree-sha1 = "a9ca650f074942cadb98d59198f93b9adc11d1ea" +uuid = "a2bd30eb-e257-5431-a919-1863eab51364" +version = "1.0.1" + +[[ImageCore]] +deps = ["ColorTypes", "Colors", "FFTW", "FixedPointNumbers", "Graphics", "MappedArrays", "OffsetArrays", "PaddedViews", "Reexport"] +git-tree-sha1 = "17a646b2b4f018d36170a99cb6a690773bb69ede" +uuid = "a09fc81d-aa75-5fe9-8630-4744c3626534" +version = "0.8.7" + +[[ImageDistances]] +deps = ["ColorVectorSpace", "Colors", "Distances", "FixedPointNumbers", "ImageCore"] +git-tree-sha1 = "aa69ce81260bcb5e950a5e3b48ccca15447c6d8c" +uuid = "51556ac3-7006-55f5-8cb3-34580c88182d" +version = "0.2.4" + +[[ImageFiltering]] +deps = ["CatIndices", "ColorVectorSpace", "Colors", "ComputationalResources", "DataStructures", "FFTViews", "FFTW", "FixedPointNumbers", "ImageCore", "LinearAlgebra", "MappedArrays", "OffsetArrays", "Requires", "StaticArrays", "Statistics", "TiledIteration"] +git-tree-sha1 = "e7be8c871a1e38cd29d1e978f7bff3feb7079025" +uuid = "6a3955dd-da59-5b1f-98d4-e7296123deb5" +version = "0.6.5" + +[[ImageQualityIndexes]] +deps = ["ColorVectorSpace", "ImageCore", "ImageDistances", "ImageFiltering", "MappedArrays", "OffsetArrays", "Statistics"] +git-tree-sha1 = "f6a2f42210340ba3680d94241aa7401b03078842" +uuid = "2996bd0c-7a13-11e9-2da2-2f5ce47296a9" +version = "0.1.2" + +[[ImageTools]] +deps = ["AlgTools", "ColorTypes", "FileIO", "GR", "OffsetArrays", "Printf", "Setfield"] +path = "../ImageTools" +uuid = "b548cc0d-4ade-417e-bf62-0e39f9d2eee9" +version = "0.1.0" + +[[IntelOpenMP_jll]] +deps = ["Libdl", "Pkg"] +git-tree-sha1 = "fb8e1c7a5594ba56f9011310790e03b5384998d6" +uuid = "1d5cc7b8-4909-519e-a0f8-d0f5ad9712d0" +version = "2018.0.3+0" + +[[InteractiveUtils]] +deps = ["Markdown"] +uuid = "b77e0a4c-d291-57a0-90e8-8db25a27a240" + +[[IntervalSets]] +deps = ["Dates", "Statistics"] +git-tree-sha1 = "4214b48a62eb8f2c292b2ee34a508c256c0cdbc9" +uuid = "8197267c-284f-5f27-9208-e0e47529a953" +version = "0.3.2" + +[[IterTools]] +git-tree-sha1 = "05110a2ab1fc5f932622ffea2a003221f4782c18" +uuid = "c8e1da08-722c-5040-9ed9-7db0dc04731e" +version = "1.3.0" + +[[LibGit2]] +uuid = "76f85450-5226-5b5a-8eaa-529ad045b433" + +[[Libdl]] +uuid = "8f399da3-3557-5675-b5ff-fb832c97cbdb" + +[[LinearAlgebra]] +deps = ["Libdl"] +uuid = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" + +[[Logging]] +uuid = "56ddb016-857b-54e1-b83d-db4d58db5568" + +[[MKL_jll]] +deps = ["Libdl", "Pkg"] +git-tree-sha1 = "61069ae718b8ab1e325bbfb4e5268902e7ea08e3" +uuid = "856f044c-d86e-5d09-b602-aeab76dc8ba7" +version = "2019.0.117+0" + +[[MacroTools]] +deps = ["DataStructures", "Markdown", "Random"] +git-tree-sha1 = "e2fc7a55bb2224e203bbd8b59f72b91323233458" +uuid = "1914dd2f-81c6-5fcd-8719-6d5c9610ff09" +version = "0.5.3" + +[[MappedArrays]] +git-tree-sha1 = "e2a02fe7ee86a10c707ff1756ab1650b40b140bb" +uuid = "dbb5928d-eab1-5f90-85c2-b9b0edb7c900" +version = "0.2.2" + +[[Markdown]] +deps = ["Base64"] +uuid = "d6f4376e-aef5-505a-96c1-9c027394607a" + +[[Missings]] +deps = ["DataAPI"] +git-tree-sha1 = "de0a5ce9e5289f27df672ffabef4d1e5861247d5" +uuid = "e1d29d7a-bbdc-5cf2-9ac0-f12de2c33e28" +version = "0.4.3" + +[[Mmap]] +uuid = "a63ad114-7e13-5084-954f-fe012c677804" + +[[NaNMath]] +git-tree-sha1 = "928b8ca9b2791081dc71a51c55347c27c618760f" +uuid = "77ba4419-2d1f-58cd-9bb1-8ffee604a2e3" +version = "0.3.3" + +[[OffsetArrays]] +git-tree-sha1 = "87d0a91efe29352d5caaa271ae3927083c096e33" +uuid = "6fe1bfb0-de20-5000-8ca7-80f57d26f881" +version = "0.11.4" + +[[OpenSpecFun_jll]] +deps = ["Libdl", "Pkg"] +git-tree-sha1 = "65f672edebf3f4e613ddf37db9dcbd7a407e5e90" +uuid = "efe28fd5-8261-553b-a9e1-b2916fc3738e" +version = "0.5.3+1" + +[[OrderedCollections]] +deps = ["Random", "Serialization", "Test"] +git-tree-sha1 = "c4c13474d23c60d20a67b217f1d7f22a40edf8f1" +uuid = "bac558e1-5e72-5ebc-8fee-abe8a469f55d" +version = "1.1.0" + +[[PaddedViews]] +deps = ["OffsetArrays", "Test"] +git-tree-sha1 = "7da3e7e1a58cffbf10177553ae95f17b92516912" +uuid = "5432bcbf-9aad-5242-b902-cca2824c8663" +version = "0.4.2" + +[[Pkg]] +deps = ["Dates", "LibGit2", "Libdl", "Logging", "Markdown", "Printf", "REPL", "Random", "SHA", "UUIDs"] +uuid = "44cfe95a-1eb2-52ea-b672-e2afdf69b78f" + +[[Printf]] +deps = ["Unicode"] +uuid = "de0858da-6303-5e67-8744-51eddeeeb8d7" + +[[QuartzImageIO]] +deps = ["ColorTypes", "ColorVectorSpace", "FileIO", "FixedPointNumbers", "ImageCore", "Libdl"] +git-tree-sha1 = "813432448a8f06a1316943c3456dcf043afacd96" +uuid = "dca85d43-d64c-5e67-8c65-017450d5d020" +version = "0.6.0" + +[[REPL]] +deps = ["InteractiveUtils", "Markdown", "Sockets"] +uuid = "3fa0cd96-eef1-5676-8a61-b3b8758bbffb" + +[[Random]] +deps = ["Serialization"] +uuid = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c" + +[[RangeArrays]] +deps = ["Compat"] +git-tree-sha1 = "d925adfd5b01cb46fde89dc9548d167b3b136f4a" +uuid = "b3c3ace0-ae52-54e7-9d0b-2c1406fd6b9d" +version = "0.3.1" + +[[Reexport]] +deps = ["Pkg"] +git-tree-sha1 = "7b1d07f411bc8ddb7977ec7f377b97b158514fe0" +uuid = "189a3867-3050-52da-a836-e630ba90ab69" +version = "0.2.0" + +[[Requires]] +deps = ["UUIDs"] +git-tree-sha1 = "999513b7dea8ac17359ed50ae8ea089e4464e35e" +uuid = "ae029012-a4dd-5104-9daa-d747884805df" +version = "1.0.0" + +[[SHA]] +uuid = "ea8e919c-243c-51af-8825-aaa63cd721ce" + +[[Serialization]] +uuid = "9e88b42a-f829-5b0c-bbe9-9e923198166b" + +[[Setfield]] +deps = ["ConstructionBase", "Future", "MacroTools", "Requires"] +git-tree-sha1 = "ed5045722fcdacf263a90fb9eb9f258598ccebac" +uuid = "efcf1570-3423-57d1-acb7-fd33fddbac46" +version = "0.5.4" + +[[SharedArrays]] +deps = ["Distributed", "Mmap", "Random", "Serialization"] +uuid = "1a1011a3-84de-559e-8e89-a11a2f7dc383" + +[[Sockets]] +uuid = "6462fe0b-24de-5631-8697-dd941f90decc" + +[[SortingAlgorithms]] +deps = ["DataStructures", "Random", "Test"] +git-tree-sha1 = "03f5898c9959f8115e30bc7226ada7d0df554ddd" +uuid = "a2af1166-a08f-5f64-846c-94a0d3cef48c" +version = "0.3.1" + +[[SparseArrays]] +deps = ["LinearAlgebra", "Random"] +uuid = "2f01184e-e22b-5df5-ae63-d93ebab69eaf" + +[[SpecialFunctions]] +deps = ["OpenSpecFun_jll"] +git-tree-sha1 = "268052ee908b2c086cc0011f528694f02f3e2408" +uuid = "276daf66-3868-5448-9aa4-cd146d93841b" +version = "0.9.0" + +[[StaticArrays]] +deps = ["LinearAlgebra", "Random", "Statistics"] +git-tree-sha1 = "5a3bcb6233adabde68ebc97be66e95dcb787424c" +uuid = "90137ffa-7385-5640-81b9-e52037218182" +version = "0.12.1" + +[[Statistics]] +deps = ["LinearAlgebra", "SparseArrays"] +uuid = "10745b16-79ce-11e8-11f9-7d13ad32a3b2" + +[[StatsBase]] +deps = ["DataAPI", "DataStructures", "LinearAlgebra", "Missings", "Printf", "Random", "SortingAlgorithms", "SparseArrays", "Statistics"] +git-tree-sha1 = "c53e809e63fe5cf5de13632090bc3520649c9950" +uuid = "2913bbd2-ae8a-5f71-8c99-4fb6c76f3a91" +version = "0.32.0" + +[[Test]] +deps = ["Distributed", "InteractiveUtils", "Logging", "Random"] +uuid = "8dfed614-e22c-5e08-85e1-65c5234f0b40" + +[[TestImages]] +deps = ["AxisArrays", "FileIO", "Pkg", "ZipFile"] +git-tree-sha1 = "8e4a38c6fed6c24f0e56e0864207973a0b48aa83" +uuid = "5e47fb64-e119-507b-a336-dd2b206d9990" +version = "1.0.0" + +[[TiledIteration]] +deps = ["OffsetArrays", "Test"] +git-tree-sha1 = "58f6f07d3b54a363ec283a8f5fc9fb4ecebde656" +uuid = "06e1c1a7-607b-532d-9fad-de7d9aa2abac" +version = "0.2.3" + +[[UUIDs]] +deps = ["Random", "SHA"] +uuid = "cf7118a7-6976-5b1a-9a39-7adc72f591a4" + +[[Unicode]] +uuid = "4ec0a83e-493e-50e2-b9ac-8f72acf5a8f5" + +[[ZipFile]] +deps = ["BinaryProvider", "Libdl", "Printf"] +git-tree-sha1 = "580ce62b6c14244916cc28ad54f8a2e2886f843d" +uuid = "a5390f91-8eb1-5f08-bee0-b1d1ffed6cea" +version = "0.8.3"
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/Project.toml Tue Apr 07 14:19:48 2020 -0500 @@ -0,0 +1,16 @@ +name = "PredictPDPS" +uuid = "172afc5a-025d-11ea-0991-59a8fed80a6d" +author = ["Tuomo Valkonen <tuomov@iki.fi>"] + +[deps] +AlgTools = "c46e2e78-5339-41fd-a966-983ff60ab8e7" +ColorTypes = "3da002f7-5984-5a60-b8a6-cbb66c0b333f" +DelimitedFiles = "8bb1440f-4735-579b-a4ab-409b98df4dab" +FileIO = "5789e2e9-d7fb-5bc7-8068-2c6fae9b9549" +GR = "28b8d3ca-fb5f-59d9-8090-bfdbd6d07a71" +ImageQualityIndexes = "2996bd0c-7a13-11e9-2da2-2f5ce47296a9" +ImageTools = "b548cc0d-4ade-417e-bf62-0e39f9d2eee9" +Printf = "de0858da-6303-5e67-8744-51eddeeeb8d7" +QuartzImageIO = "dca85d43-d64c-5e67-8c65-017450d5d020" +Setfield = "efcf1570-3423-57d1-acb7-fd33fddbac46" +TestImages = "5e47fb64-e119-507b-a336-dd2b206d9990"
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/README.md Tue Apr 07 14:19:48 2020 -0500 @@ -0,0 +1,36 @@ + +# Julia codes for “Predictive online optimisation with applications to optical flow” + +These are the Julia codes for the optical flow experiments of the manuscript _“Predictive online optimisation with applications to optical flow”_ by [Tuomo Valkonen](https://tuomov.iki.fi) ([arXiv:2002.03053](https://arxiv.org/abs/2002.03053)). + +## Prerequisites + +These codes were written for Julia 1.3. The Julia package prequisites are from November 2019 when our experiments were run, and have not been updated to maintain the same environment we used to do the experiments in the manuscript. You may get Julia from [julialang.org](https://julialang.org/). + +## Using + +Navigate your unix shell to the directory containing this `README.md` and then run: + + $ julia --project=. + +The first time doing this, to ensure all the dependencies are installed, run + + $ ]instantiate + +Afterwards in the Julia shell, type: + + > using PredictPDPS + +This may take a while as Julia precompiles the code. Then, to generate all the experiments in the manuscript, run: + + > batchrun_article() + +This will save the results under `img/`. +To see the experiments running visually, and not save the results, run + + > demo_known1() + +or any of `demo_XY()`, where `X`=1,2,3 and `Y`=`known`,`unknown`. +Further parameters and experiments are available via `run_experiments`. See the source code for details. + +To run the data generation multi-threadeadly parallel to the algorithm, set the `JULIA_NUM_THREADS` environment variable to a number larger than one.
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/src/Algorithm.jl Tue Apr 07 14:19:48 2020 -0500 @@ -0,0 +1,148 @@ +#################################################################### +# 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 + +
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/src/AlgorithmBoth.jl Tue Apr 07 14:19:48 2020 -0500 @@ -0,0 +1,239 @@ +###################################################################### +# 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 + +
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/src/AlgorithmBothCumul.jl Tue Apr 07 14:19:48 2020 -0500 @@ -0,0 +1,139 @@ +###################################################################### +# 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 + +
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/src/AlgorithmBothGreedyV.jl Tue Apr 07 14:19:48 2020 -0500 @@ -0,0 +1,167 @@ +###################################################################### +# 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 + +
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/src/AlgorithmBothMulti.jl Tue Apr 07 14:19:48 2020 -0500 @@ -0,0 +1,279 @@ +###################################################################### +# Predictive online PDPS for optical flow with unknown velocity field +###################################################################### + +__precompile__() + +module AlgorithmBothMulti + +identifier = "pdps_unknownmulti" + +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 + +
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/src/AlgorithmBothMulti.jl.orig Tue Apr 07 14:19:48 2020 -0500 @@ -0,0 +1,277 @@ +###################################################################### +# 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 + +
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/src/AlgorithmBothNL.jl Tue Apr 07 14:19:48 2020 -0500 @@ -0,0 +1,242 @@ +###################################################################### +# 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 + +
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/src/AlgorithmFB.jl Tue Apr 07 14:19:48 2020 -0500 @@ -0,0 +1,141 @@ +#################################################################### +# Predictive online PDPS for optical flow with known velocity field +#################################################################### + +__precompile__() + +module AlgorithmFB + +identifier = "fb_known" + +using Printf + +using AlgTools.Util +import AlgTools.Iterate +using ImageTools.Gradient + +using ..OpticalFlow: Image, + ImageSize, + flow! + +######################### +# Iterate initialisation +######################### + +function init_rest(x::Image) + imdim=size(x) + + y = zeros(2, imdim...) + Δx = copy(x) + Δy = copy(y) + ỹ = copy(y) + y⁻ = copy(y) + + return x, y, Δx, Δy, ỹ, y⁻ +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.ρ + τ₀, τ̃₀ = params.τ₀, params.τ̃₀ + + R_K² = ∇₂_norm₂₂_est² + τ̃ = τ̃₀/R_K² + τ = τ₀ + + ###################### + # Initialise iterates + ###################### + + x, y, Δx, Δy, ỹ, y⁻ = 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 + else + # Δx is a temporary storage variable of correct dimensions + flow!(x, v_known, Δx) + end + + ################################################################## + # We need to do forward–backward step on min_x |x-b|^2/2 + α|∇x|. + # The forward step is easy, the prox requires solving the predual + # problem of a problem similar to the original. + ################################################################## + + @. x = x-τ*(x-b) + + ############## + # Inner FISTA + ############## + + t = 0 + # Move step length from proximal quadratic term into L1 term. + α̃ = α*τ + @. ỹ = y + for i=1:params.fb_inner_iterations + ∇₂ᵀ!(Δx, ỹ) + @. Δx .-= x + ∇₂!(Δy, Δx) + @. y⁻ = y + @. y = (ỹ - τ̃*Δy)/(1 + τ̃*ρ/α̃) + proj_norm₂₁ball!(y, α̃) + t⁺ = (1+√(1+4*t^2))/2 + @. ỹ = y+((t-1)/t⁺)*(y-y⁻) + t = t⁺ + end + + ∇₂ᵀ!(Δx, y) + @. x = x - Δx + + ################################ + # Give function value if needed + ################################ + v = verbose() do + ∇₂!(Δy, x) + value = norm₂²(b-x)/2 + α*γ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/AlgorithmFBDual.jl Tue Apr 07 14:19:48 2020 -0500 @@ -0,0 +1,112 @@ +#################################################################### +# Predictive online PDPS for optical flow with known velocity field +#################################################################### + +__precompile__() + +module AlgorithmFBDual + +using Printf + +using AlgTools.Util +import AlgTools.Iterate +using ImageTools.Gradient + +using ..OpticalFlow: Image, + ImageSize, + flow! + +######################### +# Iterate initialisation +######################### + +function init_rest(x::Image) + imdim=size(x) + + y = zeros(2, imdim...) + Δx = copy(x) + Δy = copy(y) + + return x, y, Δx, Δy +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.ρ + τ₀, τ̃₀ = params.τ₀, params.τ̃₀ + + R_K² = ∇₂_norm₂₂_est² + τ = τ₀/R_K² + + ###################### + # Initialise iterates + ###################### + + x, y, Δx, Δy = init_iterates(dim) + + #################### + # Run the algorithm + #################### + + v = iterate(params) do verbose :: Function, + b :: Image, + v_known :: DisplacementT, + 🚫unused_b_next :: Image + + ################## + # Prediction step + ################## + + # Δx is a temporary storage variable of correct dimensions + flow!(@view(y[1,:, :]), Δx, v_known) + flow!(@view(y[2,:, :]), Δx, v_known) + + ∇₂ᵀ!(Δx, y) + @. x = b - Δx + ∇₂!(Δy, x) + @. y = (y - τ*Δy)/(1 + τ*ρ/α) + proj_norm₂₁ball!(y, α) + + ################################ + # Give function value if needed + ################################ + v = verbose() do + ∇₂ᵀ!(Δx, y) + @. x = b - Δx + ∇₂!(Δy, x) + value = norm₂²(b-x)/2 + α*γnorm₂₁(Δy, ρ) + value, x, [NaN, NaN], nothing + end + + v + end + + @warn "Using old x value. Better approach unimplemented as this algorithm is not good." + # ∇₂ᵀ!(Δx, y) + # @. x = b - Δx + + return x, y, v +end + +end # Module + +
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/src/ImGenerate.jl Tue Apr 07 14:19:48 2020 -0500 @@ -0,0 +1,169 @@ +################### +# 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 + +############## +# Our exports +############## + +export ImGen, + OnlineData, + imgen_square, + imgen_shake + +################## +# 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 + +################### +# 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(2) + elseif params.shaketype == :disk + return () -> begin + θ = 2π*rand(Float64) + r = params.shake*√(rand(Float64)) + return [r*cos(θ), r*sin(θ)] + end + elseif params.shaketype == :circle + return () -> begin + θ = 2π*rand(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) + +################ +# 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(sz...) + v = v_true.*(1.0 .+ params.shake_noise_level.*randn(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) + + nextv = shake(params) + v_true = nextv() + v_cumul = copy(v_true) + + while true + # Extract subwindow of original image and add noise + b_true = zeros(sz...) + extract_subimage!(b_true, im, v_cumul; threads=true) + b = b_true .+ params.noise_level.*randn(sz...) + v = v_true.*(1.0 .+ params.shake_noise_level.*randn(size(v_true)...)) + # Pass data to iteration routine + data = OnlineData{DisplacementConstant}(b_true, b, v, v_true, v_cumul) + if !put_unless_closed!(datachannel, data) + return + end + # Next step shake + v_true = nextv() + v_cumul .+= v_true + end +end + +function imgen_shake(imname, sz) + im = Float64.(Gray.(TestImages.testimage(imname))) + dynrange = maximum(im) + return ImGen(curry(generate_shake_image, im, sz), sz, 1, dynrange, + "$(imname)$(sz[1])x$(sz[2])") +end + +end # Module
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/src/OpticalFlow.jl Tue Apr 07 14:19:48 2020 -0500 @@ -0,0 +1,280 @@ +################################ +# 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=true) + end + + # We already filtered b in the previous step (b_next in that step) + b_filt = b_next_filt==nothing ? f(b) : b_next_filt + b_next_filt = f(b_next) + + return b_filt, b_next_filt +end + +end # Module
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/src/OpticalFlow.jl.orig Tue Apr 07 14:19:48 2020 -0500 @@ -0,0 +1,280 @@ +################################ +# 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
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/src/PredictPDPS.jl Tue Apr 07 14:19:48 2020 -0500 @@ -0,0 +1,508 @@ +################## +# Our main module +################## + +__precompile__() + +module PredictPDPS + +######################## +# Load external modules +######################## + +using Printf +using FileIO +#using JLD2 +using Setfield +using ImageQualityIndexes: psnr, ssim +using DelimitedFiles +import GR + +using AlgTools.Util +using AlgTools.StructTools +using AlgTools.LinkedLists +using AlgTools.Comms +using ImageTools.Visualise: secs_ns, grayimg, do_visualise +using ImageTools.ImFilter: gaussian + +##################### +# Load local modules +##################### + +include("OpticalFlow.jl") +include("ImGenerate.jl") +include("Algorithm.jl") +include("AlgorithmBoth.jl") +include("AlgorithmBothGreedyV.jl") +include("AlgorithmBothCumul.jl") +include("AlgorithmBothMulti.jl") +include("AlgorithmBothNL.jl") +include("AlgorithmFB.jl") +include("AlgorithmFBDual.jl") + +import .Algorithm, + .AlgorithmBoth, + .AlgorithmBothGreedyV, + .AlgorithmBothCumul, + .AlgorithmBothMulti, + .AlgorithmBothNL, + .AlgorithmFB, + .AlgorithmFBDual + +using .ImGenerate +using .OpticalFlow: DisplacementFull, DisplacementConstant + +############## +# Our exports +############## + +export run_experiments, + batchrun_article, + demo_known1, + demo_known2, + demo_known3, + demo_unknown1, + demo_unknown2, + demo_unknown3 + +################################### +# Parameterisation and experiments +################################### + +struct Experiment + mod :: Module + DisplacementT :: Type + imgen :: ImGen + params :: NamedTuple +end + +function Base.show(io::IO, e::Experiment) + displacementname(::Type{DisplacementFull}) = "DisplacementFull" + displacementname(::Type{DisplacementConstant}) = "DisplacementConstant" + print(io, " + mod: $(e.mod) + DisplacementT: $(displacementname(e.DisplacementT)) + imgen: $(e.imgen.name) $(e.imgen.dim[1])×$(e.imgen.dim[2]) + params: $(e.params ⬿ (kernel = "(not shown)",)) + ") +end + +const default_save_prefix="img/" + +const default_params = ( + ρ = 0, + verbose_iter = 100, + maxiter = 10000, + save_results = true, + save_images = true, + save_images_iters = Set([1, 2, 3, 5, + 10, 25, 30, 50, + 100, 250, 300, 500, + 1000, 2500, 3000, 5000, + 10000]), + pixelwise_displacement=false, + dual_flow = true, + prox_predict = true, + handle_interrupt = true, + init = :zero, + plot_movement = false, +) + +const square = imgen_square((200, 300)) +const lighthouse = imgen_shake("lighthouse", (200, 300)) + +const p_known₀ = ( + noise_level = 0.5, + shake_noise_level = 0.05, + shake = 2, + α = 1, + ρ̃₀ = 1, + σ̃₀ = 1, + δ = 0.9, + σ₀ = 1, + τ₀ = 0.01, +) + +const p_unknown₀ = ( + noise_level = 0.3, + shake_noise_level = 0.05, + shake = 2, + α = 0.2, + ρ̃₀ = 1, + σ̃₀ = 1, + σ₀ = 1, + δ = 0.9, + λ = 1, + θ = (300*200)*100^3, + kernel = gaussian((3, 3), (11, 11)), + timestep = 0.5, + displacement_count = 100, + τ₀ = 0.01, +) + +const experiments_pdps_known = ( + Experiment(Algorithm, DisplacementConstant, lighthouse, + p_known₀ ⬿ (phantom_ρ = 0,)), + Experiment(Algorithm, DisplacementConstant, lighthouse, + p_known₀ ⬿ (phantom_ρ = 100,)), + Experiment(Algorithm, DisplacementConstant, square, + p_known₀ ⬿ (phantom_ρ = 0,)) +) + +const experiments_pdps_unknown_multi = ( + Experiment(AlgorithmBothMulti, DisplacementConstant, lighthouse, + p_unknown₀ ⬿ (phantom_ρ = 0,)), + Experiment(AlgorithmBothMulti, DisplacementConstant, lighthouse, + p_unknown₀ ⬿ (phantom_ρ = 100,)), + Experiment(AlgorithmBothMulti, DisplacementConstant, square, + p_unknown₀ ⬿ (phantom_ρ = 0,)), +) + +const experiments_fb_known = ( + Experiment(AlgorithmFB, DisplacementConstant, lighthouse, + p_known₀ ⬿ (τ̃₀=0.9, fb_inner_iterations = 10)), +) + +const experiments_all = Iterators.flatten(( + experiments_pdps_known, + experiments_pdps_unknown_multi, + experiments_fb_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))" +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 +####################### + +function demo(experiment; kwargs...) + run_experiments(;experiments=(experiment,), + save_results=false, + save_images=false, + visualise=true, + recalculate=true, + verbose_iter=10, + fullscreen=true, + kwargs...) +end + +demo_known1 = () -> demo(experiments_pdps_known[3]) +demo_known2 = () -> demo(experiments_pdps_known[1]) +demo_known3 = () -> demo(experiments_pdps_known[2]) +demo_unknown1 = () -> demo(experiments_pdps_unknown_multi[3], plot_movement=true) +demo_unknown2 = () -> demo(experiments_pdps_unknown_multi[1], plot_movement=true) +demo_unknown3 = () -> demo(experiments_pdps_unknown_multi[2], plot_movement=true) + +function batchrun_article(kwargs...) + run_experiments(;experiments=experiments_all, + save_results=true, + save_images=true, + visualise=false, + recalculate=false, + 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], + psnr(x, d.b_true), + ssim(x, d.b_true), + psnr(d.b_noisy, d.b_true), + 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 + +############### +# Precompiling +############### + +# precompile(Tuple{typeof(GR.drawimage), Float64, Float64, Float64, Float64, Int64, Int64, Array{UInt32, 2}}) +# precompile(Tuple{Type{Plots.Plot{T} where T<:RecipesBase.AbstractBackend}, Plots.GRBackend, Int64, Base.Dict{Symbol, Any}, Base.Dict{Symbol, Any}, Array{Plots.Series, 1}, Nothing, Array{Plots.Subplot{T} where T<:RecipesBase.AbstractBackend, 1}, Base.Dict{Any, Plots.Subplot{T} where T<:RecipesBase.AbstractBackend}, Plots.EmptyLayout, Array{Plots.Subplot{T} where T<:RecipesBase.AbstractBackend, 1}, Bool}) +# precompile(Tuple{typeof(Plots._plot!), Plots.Plot{Plots.GRBackend}, Base.Dict{Symbol, Any}, Tuple{Array{ColorTypes.Gray{Float64}, 2}}}) + +end # Module
